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ABSTRACT 


A  study  was  made  of  the  feasibility  of  an  uncooled  throat  for  a 
large  hypersonic  wind  tunnel  facility  using  currently  available  mater¬ 
ials.  Maximum  fullscale  stagnation  conditions  would  be:  2000  psi, 
4400°R,  1500  lb/ sec  air  flow,  and  throat  diameter  10.5-inches.  The 
basic  throat  concept  was  that  of  a  ceramic  insulation  layer,  composed 
of  small  pieces,  that  would  form  a  protective  liner  within  a  metal 
structure.  High  resistance  to  thermal  spalling  was  the  material  char¬ 
acteristic  of  greatest  importance.  Tests  were  made  of  several  zirconia 
materials  and  two  zirconium  diboride  compositions  by  exposing  them  to 
hot  air  flow  in  a  sonic  throat  at  maximum  conditions  of  800  psi  and 
3550°R.  Behavior  of  the  zirconia  materials  ranged  from  minor  crack¬ 
ing  to  complete  fragmentation.  The  zirconium-diborides  did  not  crack 
and  were  oxidation  resistant  at  these  conditions.  In  addition,  the  ther¬ 
mal  stress  distribution  was  studied  for  the  individual  blocks  that  would 
form  the  throat  insulation.  For  this  purpose  the  three-dimensional 
stress  distribution  was  calculated  for  mechanically  unrestrained  blocks 
having  one -dimensional  temperature  distributions.  Effects  of  tempera¬ 
ture  distribution,  block  size  and  block  shape  were  determined.  The 
computer  program  is  included  with  the  report.  It  was  concluded  that 
currently  available  materials  are  not  satisfactory  for  a  throat  that 
would  be  used  with  no  cooling.  Some  of  the  materials  tested  may  be 
satisfactory  if  the  thermal  shock  conditions  were  reduced  by  use  of 
film  cooling  (less  than  that  required  for  a  back-side  cooled  throat)  and 
by  preheating  with  a  flow  of  air  through  the  throat  during  heater  pres¬ 
surization.  A  major  problem  in  use  of  zirconia  would  be  the  attachment 
of  the  insulation  layer  to  the  backup  structure.  Use  of  zirconium- 
diboride  would  require  a  design  concept  that  would  be  compatible  with 
its  high  thermal  conductivity.  Both  materials  might  be  used  to  advan¬ 
tage  in  a  single  design. 
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SECTION  I 
INTRODUCTION 


This  program  was  carried  out  as  part  of  the  development  effort 
for  a  large  hypersonic  wind  tunnel  facility.  This  facility  would  provide 
for  large  scale  hypersonic  aerodynamic  and  propulsion  system  testing 
under  conditions  duplicating  actual  flight.  The  necessary  high  stag¬ 
nation  temperatures  would  be  achieved  by  heating  high  pressure  air  In 
a  regenerative  heater,  using  refractory  materials  for  the  matrix. 

One  of  the  nozzles  for  this  facility  would  yield  Mach  7  flow 
and  would  have  an  exit  diameter  of  12  ft.  The  corresponding  flow  rate 
would  be  1500  lb/sec  at  stagnation  conditions  of  2000  psl  and  440° F. 
The  heat  flux  (cold  wall)  at  the  throat  would  be  about  3000  BTU/ft^sec. 

The  throat  section  for  this  nozzle  would  be  one  of  the  more 
critical  components  of  the  facility  because  of  the  combination  of  high* 
thermal  loads,  high  mechanical  loads  and  large  size.  One  design  con¬ 
cept  would  be  an  all  metal  structure  with  backside  water  cooling  and  air 
film  cooling.  Preliminary  design  work  has  Indicated  that  the  air  film 
cooling  requirement  would  be  equal  to  about  25%  of  the  mainstream  flow. 
This  large  film  cooling  rate  could  significantly  reduce  the  useful  test 
region  size  and,  therefore,  alternate  solutions  are  sought. 

One  alternative  would  be  to  use  refractory  ceramics  and/or 
metals  in  a  throat  section  so  designed  that  air  film  cooling  would  not 
be  needed.  It  was  the  purpose  of  the  work  reported  herein  to  examine 
the  feasibility  of  this  uncooled  throat  concept.  Further,  the  purpose 
was  to  establish  the  feasibility  of  existing  materials ,  as  opposed  to 
the  development  of  new  materials  • 
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SECTION  II 

PRELIMINARY  CONSIDERATIONS 


2 . 1  MACH  7  NOZZLE  AND  VALVE 

At  the  time  this  work  began,  a  preliminary  design  of  the  wind 
tunnel  facility  was  available .  The  design  incorporated  a  plug  valve 
with  the  Mach  7  nozzle.  This  valve  seated  on  the  subsonic  contraction 
and  was  to  open  by  being  retracted  into  the  stilling  chamber.  The  valve 
would  allow  a  "fast  start"  of  flow  through  the  nozzle  by  remaining  closed 
until  the  heater  had  been  pressurized  to  the  desired  tunnel  stagnation 
pressure.  Alternately,  the  valve  could  remain  open  during  heater  pres¬ 
surization.  In  this  case  air  would  flow  through  the  nozzle  as  soon  as 
pressurization  began. 

Some  pertinent  dimensions  of  the  nozzle  are:  overall  length  85  ft, 
stilling  chamber  diameter  8  ft,  exit  diameter  12  ft,  throat  diameter  10  .5- 
lnches.  The  plug  valve  would  seat  1-1/2  ft  upstream  of  the  throat,  where 
the  diameter  is  about  16-inches. 

The  nozzle  was  designed  in  sections .  The  sections  of  interest 
here  are  the  three  that  formed  the  subsonic  contraction  and  throat.  These 
were:  1)  the  contraction  section,  being  the  upstream  portion  of  the  sub¬ 
sonic  contraction;  2)  the  plug  valve  seat  section,  being  a  continuation  of 
the  subsonic  contraction;  and  3)  the  throat  section.  The  throat  section 
was  8-1/2  ft  long,  with  an  upstream  diameter  of  14-inches  and  a  down¬ 
stream  diameter  of  3  4 -inches . 

Cooling  provisions  were  as  follows:  contraction  section  backside 
water  cooled;  plug  valve  seat  air  film  cooled  with  slot  located  Immediately 
upstream  of  seat;  throat  section  backside  water  and  air  film  cooled.  Mr 
for  film  cooling  of  the  plug  valve  seat  would  aid  in  cooling  the  throat 
section,  and  additional  air  would  be  introduced  through  a  slot  located 
downstream  of  the  plug  valve  seat.  The  film  cooling  requirements  at 
design  condition  (2000  psi,  4400 °R,  1500  lb/sec)  were  9%  of  the  main¬ 
stream  flow  at  the  upstream  slot  and  16%  at  the  downstream  slot  (135 
lb/sec  and  240  lb/sec,  respectively). 

The  change  to  an  un cooled  throat  would  involve  at  least  the 
replacement  of  the  8-1/2  ft  long  throat  section.  If  the  plug  valve  seat 
were  retained,  its  film  cooling  would  cause  some  reduction  in  heat  trans¬ 
fer  to  the  throat.  If  the  seat  section  were  not  retained,  or  if  a  means 
of  eliminating  entirely  all  air  film  cooling  were  found,  the  heat  transfer 
rates  would  not  be  reduced .  The  maximum  rates  would  occur  during  a 
fast  start  of  flow  with  an  initially  cold  throat.  This  case  was  assumed 
for  the  theoretical  analysis .  It  was  not  possible  to  produce  such  large 
rates  in  the  test  apparatus . 
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2.2  UNCOOLED  THROAT  CONCEPT 

The  basic  concept  for  the  uncooled  throat  is  that  of  a  metal 
structure  lined  with  a  ceramic  insulation  material.  The  metal  structure 
would  accommodate  the  mechanical  loads  from  internal  pressure  and 
other  sources .  It  would  also  provide  means  of  attachment  to  mating 
nozzle  sections  and  would  provide  the  necessary  stiffness  in  bending. 

The  ceramic  liner  could  be  either  a  continuous  surface  or  could 
be  composed  of  separate  pieces.  The  use  of  separate  pieces  decreases 
thermal  stress  levels  and  permits  accommodation  of  thermal  expansion 
through  many  joints.  This  approach  is  often  used  with  brittle  materials. 
(One  configuration  of  the  Dynasoar  nose  cap  and  leading  edge  thermal 
protection  consisted  of  many  small  tiles  of  zirconia  fastened  to  the  back¬ 
up  structure.)  It  has  the  disadvantage  of  requiring  some  method  of  hold¬ 
ing  the  insulation  blocks  in  place. 

Very  few  materials  are  available  that  can  withstand  4000°F  in  an 
oxidation  atmosphere.  Examination  of  the  possibilities  resulted  in 
zirconia  being  selected  as  the  primary  candidate.  Zirconia  has  low 
thermal  conductivity  and,  therefore,  a  layer  one  inch  thick  would  provide 
sufficient  insulation.  On  this  basis  the  size  of  the  insulation  elements 
was  selected  as  nominally  a  one  inch  cube. 

2  .3  HEAT  TRANSFER  REDUCTION 

The  thermal  shock  to  the  throat  could  be  reduced  by  a  reduction 
in  the  heat  transfer  rate.  This  could  be  achieved  by  the  use  of  film 
cooling  and/or  by  preheating  the  insulation  layer. 

Preheating  methods  can  be  divided  into  two  classes:  low  rate 
where  the  heat  is  applied  slowly  and  high  rate  where  rapid  heating  is 
used .  Low  rate  heating  will  not  produce  large  temperature  differences 
within  the  ceramics,  but  will  raise  the  ceramic  temperature  level.  This 
could  be  done  by  electric  heaters  within  or  behind  the  ceramic.  Alter¬ 
nately,  heated  gases  could  be  passed  through  the  nozzle  throat.  Still 
another  possibility  is  to  provide  flow  passages  within  the  ceramic  wall 
and  parallel  to  the.  nozzle  surface  which  could  be  used  to  duct  preheat 
gases  from  their  source  to  the  heater  or  to  a  discharge  vent. 

Whatever  the  method,  low  level  preheating  would  require  a  metal 
backup  structure  capable  of  operating  at  relatively  high  temperatures . 

For  example,  reduction  of  the  heat  transfer  rate  by  a  factor  of  two  would 
require  preheating  to  over  2000 °F,  which  is  too  high  for  most  metals. 

The  second  approach  to  preheating  would  be  to  heat  the  throat 
at  a  relatively  high  rate,  from  the  air  side.  This  could  be  done  most 
conveniently  by  starting  nozzle  flow  with  the  hot  valve  open.  That  is, 
flow  would  pass  through  the  throat  during  heater  pressurization.  The 
disadvantage  would  be  a  loss  of  facility  run  time. 
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Air  film  cooling  would  reduce  the  heat  transfer  rates  and  would 
also  reduce  the  maximum  temperature  reached  by  the  insulation  layer . 

The  effect  of  using  the  hot  valve  slot  only,  with  9%  flow,  is  shown  in 
the  table  below. 

Wall  temperatures  with  film  cooling  and  heat  transfer  coefficients 
were  calculated  by  the  methods  of  References  4  and  5,  respectively. 

In  the  case  of  heat  transfer  coefficients,  the  values  predicted  in  Refer¬ 
ence  5  were  reduced  by  a  factor  of  2/3 .  This  factor  has  been  used  to 
account  for  typical  measured  differences  between  actual  rocket  nozzles 
and  simple  hot  air  situations.  See  for  example,  Reference  6. 


Distance  Downstream 

Maximum  Wall 

Initial  Heat  Flux 
BTU/ft^sec 

From  Throat 

Temperature 

with 

without 

Feet 

°F 

film  cool. 

film  cool 

0 

2500 

1900 

3200 

2 

2900 

1450 

2000 

4 

3150 

750 

960 

6 

3250 

380 

470 

8 

3350 

200 

240 

These  values  are  based  on  stagnation  conditions  of  2000  psi 
and  4400° R.  The  initial  heat  flux  values  are  based  on  a  wall  tempera¬ 
ture  of  300 °F.  It  is  apparent  that  this  level  of  film  cooling  would  sig¬ 
nificantly  reduce  the  heat  flux  and  the  temperature  level  within  the  throat. 
Without  film  cooling,  maximum  wall  temperature  at  the  throat  would  be 
on  the  order  of  4300°  R.  The  heat  flux  levels  and  corresponding  thermal 
shock  conditions  are  still  very  high,  far  beyond  conditions  normally  con¬ 
sidered  for  ceramic  materials . 


2  .4  TEST  PROGRAM  AND  ANALYSIS 

The  thermal  stress  levels  in  the  insulation  layer  would  exceed 
the  strength  of  available  materials  under  even  the  most  favorable  condi¬ 
tions  (small  individual  blocks  with  preheating  and  film  cooling).  Thus 
the  most  important  characteristic  required  of  the  material  is  resistance 
to  thermal  spalling.  Spalling  resistance  cannot  be  determined  with 
currently  available  analytical  methods.  Most  thermal  stress  resistance 
analyses  are  based  upon  prediction  of  conditions  that  will  produce  stresses 
equal  to  the  strength,  i.e.,  the  initiation  of  cracks.  Thermal  spalling, 
on  the  other  hand,  occurs  when  the  stresses  exceed  the  strength,  cracks 
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are  Initiated,  and  the  amount  of  damage  (spalling)  depends  upon  the 
crack  propagation  characteristics  of  the  material.  Analytical  methods, 
sufficient  for  design  purposes,  are  not  yet  available  for  prediction  of 
spalling  damage.  Tests  must  be  used. 

A  special  test  setup  was  designed  and  fabricated  for  testing  of 
candidate  materials.  The  test  arrangement  simulated  a  small  portion  of 
the  facility  throat  wall.  The  test  specimens  were  small  blocks  equi¬ 
valent  in  size  to  those  that  could  be  used  in  the  fullscale  throat.  The 
specimens  formed  the  two  parallel  sides  of  a  two-dimensional  slot  throat, 
and  were  exposed  to  high  temperature,  high  pressure  air  flow.  This 
equipment,  the  tests,  and  the  results  are  described  in  Section  III. 

Experience  and  qualitative  reasoning  indicated  that  small  blocks 
would  have  greater  spall  resistance  than  large  blocks.  Information  was 
needed  on  the  influence  of  size  and  shape  of  block.  Recognizing  that 
an  analysis  of  spalling  itself  could  not  be  done,  it  was  felt  that  a  stress 
analysis  would  provide  valuable  insight. 

The  small  blocks  that  would  form  the  insulation  liner  would  be 
essentially  unrestrained  in  terms  of  mechanically  applied  loads.  That  is, 
the  thermal  stresses  would  be  much  higher  than  the  mechanical  stresses. 
Examination  of  the  literature  indicated  that  this  case  of  thermal  stress 
had  not  been  published.  Therefore,  an  analysis  was  undertaken.  The 
stress  analysis  and  computer  programming  was  performed  by  the  Illinois 
Institute  of  Technology  Research  Institute.  This  work  is  described  in 
Section  IV,  and  the  report  by  IITRI  is  presented  in  Appendices  I  and  II, 
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SECTION  III 
TEST  PROGRAM 


3 . 1  FACILITY  AND  TEST  APPARAT  US 

The  material  tests  were  conducted  at  the  FluiDyne  Medicine  Lake 
Laboratory.  Ceramic  specimens  were  tested  by  exposure  to  hot  air  flow 
from  a  zirconia  storage  heater.  The  uncooled  throat  apparatus  was  con¬ 
nected  directly  to  the  discharge  outlet  of  the  heater.  There  was  no  pro¬ 
vision  for  a  fast  start  of  flow  over  the  specimens  (no  "hot  valve"  at  the 
heater  outlet) .  Thus,  hot  air  flowed  over  the  specimens  during  heater 
pressurization. 

A  schematic  of  the  heater  system  is  shown  in  Figure  1 .  Air  is 
stored  at  5000  psi  and  throttled  manually  to  the  desired  heater  pressure. 

The  air  is  heated  by  passing  it  through  a  bed  of  zirconia  cored  brick. 

The  bed  is  heated  with  an  air-oxygen-propane  burner. 

The  uncooled  throat  test  apparatus  consisted  of  four  major  com¬ 
ponents:  a  ceramic  nozzle  entrance  ring,  a  water  cooled  heater  back 
pressure  nozzle,  a  water  cooled  test  specimen  holder,  and  a  water  cooled 
retainer  (diffuser  section) .  The  assembly  of  these  components  is  shown 
in  Figures  2  and  3 . 

The  ceramic  nozzle  entrance  ring  was  utilized  to  reduce  the  heater 
outlet  passage  to  the  2.00-inch  heater  backpressure  (slot)  nozzle  dimension. 
This  ceramic  ring  was  fabricated  from  partially  stabilized  calcia  zirconia 
by  the  castable  process . 

The  most  critical  component,  from  a  safety  viewpoint,  was  the 
watercooled  throat  (Figure  4) .  The  purpose  of  this  throat  was  to  back¬ 
pressure  the  heater.  The  throat  itself  is  a  slot,  l/10-inch  wide  by 
2-inches  high,  located  immediately  upstream  of  the  ceramic  test  specimens. 
Failure  of  the  throat  would  allow  increased  flow  through  the  heater,  and 
possibly  flotation  of  the  heater  bed. 

The  water  cooled  test  specimen  holder  (Figure  5)  was  designed  to 
hold  the  test  specimens  in  place  and  allow  minor  adjustments  to  assure 
proper  alignment  and  centering.  Since  the  shape  of  the  test  specimens 
was  changed  between  the  two  test  phases,  the  holder  internal  cavity 
was  modified  prior  to  the  second  test  phase.  The  hot  air  slot  through  the 
holder  is  formed  by  the  installed  test  specimens  on  two  sides  and  by 
the  backside  water  cooled  zirconium  copper  holder  on  the  narrow  ends 
of  the  slot. 

The  retainer,  shown  in  Figures  2  and  6,  is  the  water  cooled 
diffuser  section  and  provides  spray  cooling  to  cool  the  exhausting  hot 
air.  This  retainer  fastens  directly  to  the  large  flange  containing  the 
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backpressure  throat  and  thereby  sandwiches  the  test  specimen  holder 
between  the  large  flange  and  the  retainer.  Alignment  between  the 
three  parts  is  maintained  using  dowel  pins . 

i 

Detailed  design  of  the  metal  surfaces  exposed  to  hot  air  flow 
established  a  requirement  for  use  of  all  available  high  pressure  cooling 
water,  100  gpm  at  900  psi.  Furthermore,  it  was  necessary  to  cool 
different  components  by  passing  the  water  through  them  in  series. 

The  metal  surfaces  protected  by  the  ceramic  test  specimens  were  not 
cooled.  These  design  compromises  introduced  an  operating  limit  on 
stagnation  pressure. 

Safety  of  operation  dictated  that  failure  of  the  ceramic  speci¬ 
mens  under  test  not  endanger  the  overall  facility.  Thus,  damage  to 
metal  parts  downsteam  of  the  backpressure  throat  would  be  tolerable, 
but  loss  of  the  throat  itself  would  not  be  tolerable,  because  throat 
failure  could  cause  flotation  of  the  heater  bed. 

The  integrity  of  the  backpressure  throat  could  be  preserved, 
even  with  failure  of  the  test  specimens ,  if  the  throat  cooling  water  was 
not  affected.  However,  the  throat  cooling  water  and  the  ceramic  speci¬ 
men  holder  cooling  water  are  in  series .  Failure  of  the  ceramic  specimens 
could  result  in  overheating  of  the  specimen  holder.  More  specifically, 
O-rlngs  that  seal  the  water  passages  could  overheat  and  leak.  At 
stagnation  pressures  below  about  1000  psi  the  leakage  would  be  from 
water  to  air  and  therefore  not  harmful.  At  air  stagnation  pressures 
above  1000  psi  the  leakage  would  be  air  into  the  water  passage.  This 
would  tend  to  restrict  the  water  flow  through  both  the  ceramic  holder 
and  the  backpressure  throat. 

Thus,  1000  psi  was  the  maximum  heater  pressure  for  which  a 
failure  of  the  ceramic  specimens  could  be  sustained  without  a  substan¬ 
tial  risk  to  the  heater.  On  this  basis  an  upper  test  pressure  limit  of 
800  psi  was  established. 


3.2  TEST  SPECIMENS 

The  test  program  was  carried  out  in  two  phases,  with  some 
difference  in  the  test  specimen  shapes .  The  arrangement  used  for 
Phase  I  is  shown  in  the  sketch  below  and  in  Figure  8.  Two  assemblies 
of  nine  pieces  each  form  the  sides  of  a  slot  l/10-inch  wide  by  2-inches 
high.  The  slot  was  vertical  to  minimize  the  possibility  of  blockage  of 
flow  by  spalled  fragments.  Separate  pieces  were  used,  with  slight 
spacing,  to  accommodate  thermal  expansion  and  reduce  stresses.  Spacing 
was. provided  by  metal  shims  near  the  cold  side.  The  assemblies  were 
placed  in  the  specimen  holder  (Figure  5)  and  were  partially  held  in  place 
with  copper  dowels  (Figure  8) . 
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A  simpler  arrangement  was  used  for  the  Phase  II  tests.  The 
assemblies  were  replaced  by  one  piece  specimens  (see  Figure  10)  and 
the  dowels  were  eliminated.  The  overall  shape  was  changed  to 
trapezoidal.  Specimens  were  held  in  place  by  shim  plates  and  set 
screws  acting  on  the  slanting  surfaces  (see  Figure  5) . 
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Only  three  materials  were  tested  in  Phase  II,  yttria-zirconia 
and  two  zirconium  diboride  alloys .  Thermal  expansion  slots  were  cut 
in  the  zirconla  blocks  (.010  wide  by  l/2-inch  deep),  but  not  in  the 
diboride  blocks  (Figure  10) . 

The  materials  tested  were  available  without  need  for  special 
development.  They  are  not  all  available  on  a  large  scale.  This  is 
especially  true  of  the  wire  reinforced  zirconia  and  the  zirconium  di¬ 
boride  materials.  A  brief  description  of  each  material  is  given  In  Table 
I. 


3.3  TEST  PROCEDURES 

The  basic  procedure  was  to  start  testing  at  low  pressures  and 
temperatures,  examine  the  specimens  after  test,  and  proceed  with  more 
severe  conditions  if  warranted.  Stagnation  pressures  and  temperature 
were  recorded  during  the  entire  test  sequence  of  heater  pressurization, 
steady  state  run  and  depressurization.  1 

Stagnation  temperature  was  measured  with  a  bare  wire  Pt-6RH/ 
Pt-30  Eh  thermocouple  in  the  heater  air  outlet .  This  outlet  is  a  zirconia 
insulated  passage  having  a  3-inch  internal  diameter  and  is  2-l/2~feet 
long.  The  passage  was  preheated  before  each  run  to  reduce  temperature 
losses  and  reduce  thermal  stresses  in  the  zirconia  insulation. 

It  was  not  possible  to  observe  the  test  specimens  during  the 
runs.  They  were  located  too  far  upstream  in  the  apparatus.  It  was  pos¬ 
sible  to  observe  the  air  discharge  from  the  test  apparatus.  However, 
this  viewing  was  partially  obscured  by  Jets  of  water,  which  discharged 
into  the  airstream  for  cooling  purposes .  The  combination  of  the  partially 
obscured  viewing  and  sound  was  helpful  in  indicating  failure  of  specimens 
on  several  occasions . 

Pressurization  rates  were  slower  than  normal  because  of  the 
desire  to  detect  some  indication  of  failure  during  pressurization.  De¬ 
pressurization  times  varied  widely  because  the  heater  is  partly  depres¬ 
surized  by  opening  a  valve  in  the  heater  exhaust  stack,  allowing  air 
to  flow  out  the  bottom  of  the  heater.  This  was  not  done  on  a  regular 
schedule,  which  caused  variations. 

The  test  specimen  holder,  with  the  specimens,  was  removed 
after  each  run.  The  specimens  were  examined  without  removing  them 
from  the  holder  unless  a  change  was  needed  for  the  subsequent  test. 
During  this  between  run  period  the  heater  bed  was  reheated.  The 
specimen  holder,  with  specimens,  would  be  replaced  after  the  desired 
heater  bed  temperature  was  reached. 
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The  heater  outlet  passage  was  then  preheated  by  flowing  com¬ 
bustion  gases.  These  gases  were  withdrawn  just  upstream  of  the 
heater  backpressure  throat,  through  a  side  opening  in  the  pipe.  During 
this  period  the  test  specimens  were  cooled  by  flowing  air  over  them 
from  the  open  downstream  end. 

After  completion  of  preheating,  the  specimen  cooling  air  line 
was  removed  and  the  heater  pressurization  was  started.  If  unusual 
sights  or  sounds  were  detected  during  pressurization,  the  run  was 
terminated.  Observation  continued  during  the  run  on  the  same  basis. 


3.4  TEST  RESULTS 

All  the  tests  reported  herein  were  made  at  throat  heat  transfer 
rates  substantially  less  than  that  corresponding  to  the  fullscale  wind 
tunnel  conditions .  The  fullscale  cold  wall  heat  flux  would  be  about 
3000  BTU/ft2sec  at  stagnation  conditions  of  2000  psi  and  4400° R. 

The  cold  wall  heat  flux  values  for  the  test  varied  from  400  to  1800  BTU/ 
ft2sec. 


Cold  wall  values  apply  only  if  flow  is  established  very  quickly, 
and  with  the  wall  initially  cold,  i.e.,  room  temperature.  This  is  pos¬ 
sible  with  a  hot  valve,  as  discussed  in  Section  2.1.  Flow  through  the 
throat  during  heater  pressurization  reduces  the  heat  transfer  rates  from 
the  calculated  cold  wall  values.  In  the.following  discussion,  cold  wall 
heat  flux  values  are  cited  for  comparison  purposes.  However,  the 
actual  heat  flux  rates  during  the  tests  were  less  than  these  values,  and 
this  fact  must  be  considered  in  examining  the  results  for  application  to 
fullscale  conditions  • 

The  specimens  that  were  tested  are  identified  in  Table  I.  A 
summary  of  the  test  conditions  is  given  in  Table  II.  Photographs  of 
the  test  specimens  before  and  after  tests  are  presented  in  Figures  8 
through  17 .  The  damage  in  some  instances  is  not  apparent  from  the 
photographs,  but  is  described  below. 


3*4.1  Phase  I  Tests 

Tuncsten-Rhenium  Wire  Reinforced  Zirconia 

Specimens  1  and  2.  Runs  1,  2,  3.  Figures  8  and  11.  Cold 
wall  heat  flux  400  to  500  BTU/ft2sec. 

Oxidation  of  the  wire  reinforcement  is  a  basic  limitation  of  this 
material.  It  was  hoped  that  specimens  of  low  porosity  could  be  tested. 
These  would  have  wires  exposed  only  at  the  surface  and,  therefore, 
would  provide  the  least  exposure  to  oxidation. 
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However,  the  specimens  that  were  available  and  tested  had 
large  lamination-type  pores  that  exposed  wires  well  below  the  surface. 
Under  test  the  specimens  swelled  and  began  to  flake.  Additional 
testing  would  have  caused  spalling  and  further  oxidation. 

The  zirconia  material  was  a  low  thermal  expansion  composition, 
and  was  chosen  for  this  reason.  Low  thermal  expansion  would  reduce 
thermal  stresses  and  enhance  the  probability  of  success.  These  low 
expansion  compositions  are  discussed  further  below. 

Medium  Grain  Size  Zirconia 

Specimens  3  and  4.  Run  4.  Figures  8  and  11.  Cold  wall  heat 
flux  650  BTU/ft^sec. 

Specimen  3  was  castable  calcla-zirconia  and  Specimen  4  was 
castable  yttria -zirconia  with  maximum  grain  size  of  60  mesh.  Both 
spalled  severely  to  a  depth  of  1/8  to  l/4-inch.  The  yttria -zirconia 
is  at  the  left  side  of  the  photo  in  Figure  11. 

Fine  Grain  Dense  Zirconia 

Specimens  5  and  6.  Runs  5,  '6,  7,  11.  Figures  8  and  11.  Cold 
wall  heat  flux  450  to  1000  BTU/ft^sec. 

This  material  is  a  very  dense  (about  5%  porosity)  partially  stabil¬ 
ized  zirconia  formulated  to  have  a  low  thermal  expansion  coefficient. 

The  principal  stabilizers  are  magnesia  and  calcia.  Thus  the  compositional 
stability  at  high  temperatures  is  limited.  Whether  the  life  would  be  ade¬ 
quate  for  an  uncooled  throat  is  still  uncertain. 

This  material  is  currently  used  in  several  commercial  applications . 
It  was  selected  in  order  to  try  a  dense  material,  recognizing  that  cracks 
could  not  be  avoided .  The  low  thermal  expansion  characteristics  would 
enhance  the  probability  of  success. 

The  tests  produced  a  mosaic  pattern  of  hairline  cracks,  spaced 
about  1 /4-inch.  There  were  some  open  cracks  in  the  side  pieces  and 
spalling  of  one  downstream  piece.  Some  of  the  cracks  entered  the  sur¬ 
face  at  shallow  angles,  indicating  that  more  severe  tests  would  cause 
spalling. 

This  material  may  be  suitable  for  small  throats,  especially  if 
an  initial  heat  shock  can  be  used  to  precrack  the  material  (Reference  1) . 

In  a  large  throat,  it  is  unlikely  that  a  crack  pattern  could  be  produced 
that  would  relieve  stresses  and  not  allow  spalling. 
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Coarse  Grain  Calcla-Zlrconla 

Specimens  7 ’and  8.  Runs  8  and  15.  Figures  9,  11,  13.  Cold 
wall  heat  flux  650  and  1100  BTU/ftzsec. 

This  material  is  a  partially  stabilized  castable  zirconia.  It  is 
used  very  successfully  in  wind  tunnel  air  heaters  as  dense  insulation. 

In  this  application  it  has  shown  excellent  resistance  tt>  thermal  spalling. 
The  large  grain  size  (to  8  mesh)  is  a  disadvantage  for  the  present  appli¬ 
cation.  It  does  not  provide  a  smooth  surface,  and  it  is  difficult  to  main¬ 
tain  comers  and  edges  even  without  severe  thermal  stress  conditions. 

For  this  reason  a  finer  grain  60  mesh  material  was  also  included  in  the 
tests.  However,  its  performance  was  unsatisfactory,  (Specimen  3) . 

The  test  of  Run  8  produced  some  spalling  of  Specimen  7,  as 
shown  in  Figure  11,  but  no  spalling  of  Specimen  8.  This  indicates  that 
the  material  may  not  have  the  uniformity  needed  for  a  large  throat 
application.  A  second  test  was  made  with  Specimen  8.  This  test, 

Run  15,  did  not  produce  useful  results  because  of  partial  melting  of 
the  copper  dowels  (Figure  13) . 

Coarse  Grain  Magnesla-Calcla -Zirconia 

Specimens  9  and  10.  Runs  9,  10,  14.  Figures  9  and  12.  Cold 
wall  heat  flux  650  to  950  BTU/ftzsec, 

This  material  was  also  a  castable,  similar  to  the  calcia -zirconia 
discussed  above,  but  with  the  magnesia  and  calcia  stabilizers  also 
mentioned  earlier.  Its  behavior  was  similar  to  the  coarse  grain  calcia- 
zirconia.  Fine  cracks  were  found  on  some  of  the  pieces,  and  there 
was  a  slight  amount  of  spalling,  about  l/32-inch  deep.  This  can  be 
seen  in  the  center  piece  in  the  upper  left  hand  photo  in  Figure  12. 

Medium  Grain  Macmesla -Calcia -Zirconia 

Specimen  11.  Runs  11  and  15.  Figures  9,  12,  13.  Cold  wall 
heat  flux  1000  to  1100  BTU/ft2sec. 

This  material  differed  from  that  discussed  immediately  above 
only  in  maximum  grain  size.,  having  an  upper  size  of  60  mesh.  The 
specimen  survived  Run  11  with  no  visible  damage.  Run  15  did  not 
produce  useful  results  because  of  melting  of  the  copper  dowels. 

These  results  are  quite  different  from  those  of  Run  4,  with  the 
60  mesh  castable  calcia -zirconia  and  yttria -zirconia .  The  difference 
could  be  attributed  to  the  lower  thermal  expansion  coefficient  of  the 
1027  type  composition.  But  there  is  not  enough  information  to  determine 
the  reason  for  the  difference  in  performance. 


12 


AEDC-TR-70-92 


Pressed  Magnesia -Zlrconla 

Specimens  12  and  13.  Runs  12  and  13.  Figures  9  and  12. 

Cold  wall  heat  flux  850  to  1100  BTU/ft2Sec. 

This  material  is  another  composition  having  low  thermal  ex¬ 
pansion  characteristics.  It  was  fabricated  by  pressing,  as  was  the 
fine  grain  dense  zlrconla.  The  tests  caused  formation  of  fine  cracks, 
but  no  spalling  occurred.  It  Is  uncertain  whether  the  cracks  would 
lead  to  spalling  after  repeated  thermal  stress  cycles. 

3.4.2  Phase  II  Tests 

Yttria -Zlrconla 

Specimens  14,  15,  16  and  17.  Runs  1  and  5.  Figures  10,  14, 

16.  Cold  wall  heat  flux  1000  BTU/ft^sec. 

These  specimens  were  prepared  in  an  attempt  to  produce  a  material 
that  would  have  high  resistance  to  thermal  spalling.  A  combination  of 
fine  and  coarse  grains  were  used  in  an  effort  to  impart  resistance  to 
crack  propagation.  The  composition  was  fully  stabilized  yttria -zirconia. 

The  results  clearly  demonstrated  that  the  objective  was  not 
achieved.  In  fact,  this  material  suffered  the  most  complete  destruction 
of  any  materials  tested.  The  test  specimens  were  completely  shattered 
during  Run  1.  A  repeat  run  was  made,  which  confirmed  this  result. 

The  extent  of  the  damage  is  evident  from  the  photographs  (Figures 
14  and  16) .  An  examination  of  the  material,  after  test,  was  made  to 
determine  if  phase  inversion  had  contributed  to  the  failure.  The  mono¬ 
clinic  content  was  measured  at  three  locations,  hot  face,  center  and 
cold  face.  All  values  were  one  percent.  Hence,  the  failure  was  caused 
by  thermal  stresses,  not  by  phase  inversion. 

Zlrconlum-Dlborlde 

Two  zirconium-diboride  compositions  were  tested.  These  compo¬ 
sitions  were  selected  from  a  family  of  materials  that  have  been  developed 
by  ManLabs  under  Air  Force  Contract  AF33 (615) -3671 . 

The  materials  were  tested  in  an  effort  to  determine  whether  their 
oxidation  resistance  would  be  sufficient  for  this  application.  Oxidation 
protection  is  provided  by  a  layer  of  zirconia  which  forms  upon  exposure 
to  oxygen,  and  which  provides  a  barrier  to  further  oxidation.  If  this 
coating  does  not  spall  under  test,  thereby  exposing  parent  material  to 
catastrophic  oxidation,  the  material  could  probably  be  used  in  throat  . 
construction. 
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The  diborides  tested  have  much  higher  thermal  conductivities 
than  zirconia.  Their  conductivities  lie  in  the  range  35  to  65  BTU/ft.hr°R, 
and  therefore,  are  roughly  twice  as  high  as  that  of  carbon  steel.  The 
high  conductivity  imparts  thermal  shock  resistance,  but  does  not  pro¬ 
vide  thermal  insulation  qualitites.  The  question  of  how  this  material 
could  be  used  in  an  "uncooled"  throat  was  not  broached.  It  was  Judged 
of  first  priority  to  determine  if  the  material  itself  was  useful. 

ManLabs  material  V,  Specimens  18  and  19,  were  tested  in  Runs 
2,  3  and  9  at  cold  wall  heat  fluxes  from  850  to  1600  BTU/ft^sec.  Man- 
Labs  material  VIII,  Specimens  22  and  23,  were  tested  in  Runs  4,  6,  7 
and  8,  at  cold  wall  heat  fluxes  from  1000  to  1800  BTU/ft2sec.  The 
specimens  after  test  are  shown  in  Figures  15  and  17  • 

The  high  thermal  conductivity  of  the  test  specimens  required 
that  the  copper  specimen  holder  be  protected  by  insulation.  Layers  of 
zirconia,  l/4-inch  thick,  were  placed  between  the  test  pieces  and 
the  copper  surfaces.  In  addition,  stainless  steel  metal  shim  plates 
were  placed  outside  the  zirconia  insulators  (see  Figures  15  and  17). 
These  shim  plates  acted  as  bearing  plates  for  set  screws  that  held  the 
assemblies  in  position. 

All  four  test  specimens  survived  the  early  test  without  damage. 

A  layer  of  grayish  zirconia  was  visible  on  the  wetted  flow  surface  (see 
Figure  15). 

During  the  last  run  with  each  pair  of  specimens.  Runs  8  and  9, 
the  air  discharge  from  the  test  apparatus  became  brilliant  white,  and 
the  test  in  each  case  was  terminated .  In  each  case  the  downstream 
end  of  one  of  the  specimens  had  oxidized  (Figure  17) .  These  specimens 
were  examined  by  ManLabs  in  an  effort  to  understand  the  cause  of  the 
failure.  We  also  examined  photographs  that  were  taken  after  each  run. 

The  cause  of  the  failure  was  reaction  of  iron  oxide  with  the  pro¬ 
tective  zirconia  layer.  Iron  oxide  and  zirconia  form  a  eutectic  at  2380°F, 
which  is  well  below  the  temperature  of  the  test.  Thus  the  zirconia- 
iron  oxide  layer  would  be  sloughed  off,  exposing  the  diboride  material 
to  the  high  temperature  air  stream.  Under  these  conditions  of  high 
temperature  and  high  shear,  the  protective  zirconia  layer  apparently 
cannot  form,  allowing  continuous  oxidation  and  erosion. 

The  iron  oxide  undoubtedly  came  from  the  stainless  steel  shim 
pieces  in  the  following  way.  Examination  of  the  after  run  photographs 
showed  that  the  l/4-inch  zirconia  insulators  were  eroded  near  the 
throat  exit.  This  erosion  would  allow  flow  along  the  sides  of  the 
specimens,  which  would  continue  erosion  of  the  zirconia  insulators. 
Eventually  this  bypass  flow  would  reach  the  steel  shim  plates  •  At  this 
point  the  plates  would  oxidize,  and  the  iron  oxide  would  be  carried 
into  the  stream,  and  into  the  wake  behind  the  test  specimens.  There 
it  would  combine  with  the  zirconia  layer  on  the  diboride  specimens . 
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The  flow  conditions  and  time  required  to  produce  the  failures 
depended  upon  the  structural  integrity  of  the  down  stream  zirconia 
insulators  •  Thus  the  conditions  at  which  the  failures  occurred  bear  no 
relationship  to  maximum  conditions  that  the  zirconia  diborlde  can 
tolerate.  The  failures  do  emphasise  the  importance  of  not  contamin¬ 
ating  the  oxide  layer  on  these  materials . 

3.5  CONCLUSIONS  FROM  TEST  RESULTS 

In  examining  the  results  of  these  tests,  it  must  be  remembered 
that  the  test  conditions  were  not  as  severe  as  those  corresponding  to 
the  fullscale  wind  tunnel  facility.  Hence,  the  results  can  be  described, 
but  extrapolation  to  more  severe  conditions  is  difficult.  The  results 
and  conclusions  are  summarized  as  follows: 

1.  The  tungsten  wire  reinforced  zirconia  that  was  tested  is 
not  suitable  for  this  application  because  excessive  porosity 
allows  oxidation  of  the  wires .  Whether  a  less  porous 
material  would  be  satisfactory  is  unknown. 

2 .  The  behavior  of  different  zirconia  materials  varied  greatly, 
from  only  minor  cracking  to  complete  fragmentation . 

3 .  A  very  dense  zirconia  (about  5  percent  porosity)  in  four  tests 
‘  developed  a  pattern  of  tight  cracks,  but  did  not  spall  exces¬ 
sively.  Intentional  precracking  by  thermal  shock  could  make 
this  a  useful  material  for  small  throats .  For  a  larger  throat, 
with  a  segmented  structure,  the  extensive  cracking  would 
probably  lead  to  spalling . 

4.  Some  of  the  medium  and  coarse  grain  zirconias  of  nominally 
25  percent  porosity  showed  little  or  no  spalling  and,  there¬ 
fore,  show  promise.  These  materials  were  partially  stabilized 
with  calcia  and/or  magnesia.  The  magnesia  stabilized  materials 
have  low  thermal  expansion  characteristics  which  is  highly 
desirable.  On  the  other  hand,  these  compositions  are  not 
stable  at  high  temperatures.  Therefore,  further  information 

on  both  spalling  resistance  and  life  at  high  temperatures 
would  be  needed  to  determine  their  suitability. 

5.  Zirconium  diboride  was  tested  to  stagnation  conditions  of  800 
psi  and  3550° R  without  failure.  This  material  has  potential 
for  throat  construction  if  a  design  can  be  devised  that  would 
be  compatible  with  its  high  thermal  conductivity. 

6.  On  the  basis  of  this  work  and  other  work  with  ceramics,  it 

is  the  authors '  opinion  that  no  commerically  available  material 
meets  the  requirements  for  the  most  severe  conditions  of  an 
uncooled  throat  for  hypersonic  wind  tunnel  facilities .  Some 
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reduction  in  thermal  shock  would  be  necessary,  most 
probably  by  film  cooling,  especially  during  the  initial 
flow  period  when  heat  transfer  rates  are  high.  The 
cooling  rate  could  be  reduced  after  this  initial  flow 
period.  Preheating  by  flowing  gases  through  the  throat 
during  heater  pressurization  would  also  reduce  stresses. 
This  method  would  cause  some  reduction  in  facility  run 
time . 
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SECTION  IV 

THERMAL  STRESS  ANALYSIS 


4.1  INTRODUCTION 

In  parallel  with  the  test  program,  an  analytic  study  of  thermal 
stresses  was  performed.  This  study  consisted  basically  of  development 
of  a  computer  program  and  then  use  of  the  program  to  obtain  numerical 
solutions  for  a  variety  of  cases  based  on  the  full  scale  wind  tunnel 
throat . 


The  computer  program  was  developed  under  subcontract  by  the 
Illinois  Institute  of  Technology  Research  Institute.  Dr.  M.A.  Salmon 
of  IITRI  directed  development  of  the  program  and  running  of  the  initial 
solutions .  IITRI  Final  Report  M620 1 ,  Analysis  of  Thermal  Stresses  in 
Rectangular  Blocks .  summarizes  their  work.  This  report  is  included 
as  Appendix  I . 

The  necessity  of  undertaking  this  analytical  study  stems  from  a 
basic  scarcity  of  thermal  stress  information  on  unrestrained  solids.  This 
class  of  problems  is  particularly  difficult  to  analyse  (especially  so  be¬ 
cause  the  stress  field  is  three  dimensional);  hence  almost  no  information 
is  available  in  the  literature.  Developments  in  finite -element  theory 
have  resulted  in  the  ability  to  analyse  thermal  stresses  in  externally 
unrestrained  solid  bodies  such  as  the  ceramic  elements  of  the  present 
uncooled  throat  concept. 

The  model  used  in  the  analysis  is  a  right  parallelepiped  (Figure 
21)  with  a  one  dimensional  (in  the  direction  of  one  of  the  major  axes) 
temperature  distribution.  In  the  analysis  the  parallelepiped  is  subdi¬ 
vided  into  a  larger  number  of  smaller  parallelepipeds  which  are  in  fact 
the  finite  elements .  The  material  properties  can  be  allowed  to  vary 
with  temperature. 

In  applying  this  analysis  to  high  temperature  situations,  it  is 
essential  to  keep  in  mind  that  the  analysis  is  based  strictly  on  elastic 
theory.  Effects  of  plastic  deformation  in  reducing  stress  levels,  and 
in  causing  residual  stresses  are  ignored. 

The  underlying  purpose  in  performing  the  analytical  study  was 
to  improve  the  basic  understanding  of  the  stress  fields  within  the  ceramic 
elements.  Qualitative  (and  to  a  useful  extent  quantitative)  information 
can  be  gained,  for  instance  on  the  relationships  of  stresses  to  temper¬ 
ature  profile,  geometry,  material  properties,  etc.  Knowledge  of  these 
relationships  can  be  used  to  guide  selection  of  shapes,  operating 
procedure,  etc. 
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4.2  THERMAL  STRESS  PARAMETERS 

As  described  in  other  parts  of  this  report,  the  uncooled  throat 
concept  pursued  in  this  program  is  one  of  a  large  number  of  small 
ceramic  elements  as  opposed  to  a  monolithic  throat  structure.  Each 
element  is  to  be  held  in  place  with  essentially  no  external  restraints 
to  thermal  growth  of  the  element.  Even  though  external  restraints  are 
assumed  not  to  exist,  large  thermal  stresses  can  be  produced  in  a 
solid  element  such  as  the  parallelepiped  configuration  employed. 

Thermal  stresses  in  such  a  body  arise  due  to  nonlinear  temperature 
variations  through  it.  An  elastic  body  with  constant  properties 
experiences  no  thermal  stresses  if  it  is  subjected  to  a  linear  (in 
cartesian  coordinates)  temperature  variation.*  When  the  temperature 
distribution  is  nonlinear,  however,  thermal  stresses  arise.  Reference 
3  contains  derivations  and  proofs  of  the  above  relationships . 

Broadly  speaking,  the  magnitude  of  the  thermal  stresses  depends 
on  the  amount  by  which  the  actual  temperature  distribution  deviates  from 
being  linear.  Furthermore,  the  quality  of  the  stress  field  (distribution 
of  compressive  or  tensile  stresses)  will  depend  on  whether  the  curva¬ 
ture  is  positive  or  negative  (i.e.,  the  algebraic  sign  of  the  second 
derivative  of  the  temperature  distribution)  • 

Again  in  general  terms,  the  stresses  within  a  block  will  depend 
on: 

1 .  Temperature  distribution  through  the  block  • 

2 .  Overall  size  and  shape  of  the  block  • 

3.  Material  properties;  modulus  of  elasticity,  Poissons 
ratio  and  coefficient  of  thermal  expansion. 


4.2.1  Temperature  Distributions 

.  Figure  18  shows  the  time-temperature  history  for  blocks  of 
the  facility  uncooled  throat.  The  calculation  assumes  an  instantaneous 
start  of  flow  at  2000  psi  and  4400 °R.  It,  therefore,  assumes  use  of 


*Another  useful  property  of  thermal  stress  behavior  is  that  superposition 
of  thermal  stresses  resulting  from  different  temperature  profiles  applies. 
Since  a  linear  distribution  induces  no  stresses,  any  linear  temperature 
distribution  can  be  superimposed  on  an  existing  temperature  profile  with¬ 
out  altering  the  stresses.  In  this  way  stress  fields  can  sometimes  be 
more  easily  visualized. 
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a  hot  valve  and  is  the  most  severe  thermal  shock  conditions  that  the 
facility  could  produce.  It  was  selected  for  calculation  purposes  even 
through  the  stress  levels  reach  values,  in  some  instances,  well  above 
the  strength  of  available  materials.  Stresses  that  exceed  the  material 
strength  do  not  necessarily  cause  spalling.  Furthermore,  the  primary 
purpose  of  the  analysis  was  to  gain  an  understanding  of  the  stress  field. 

Note  that  as  the  time  increases,  the  hot  wall  temperature  in¬ 
creases,  but  the  value  of  second  derivative  of  the  temperature  curve 
decreases.  Therefore,  stresses  near  the  heated  surface  will  decrease 
with  time.  After  a  run  is  completed,  the  inner  surface  of  the  throat  ■ 
will  begin  to  cool,  with  the  result  that  the  temperature  curve  near  this 
surface  will  become  reversed  in  curvature.  Stresses  near  the  surface 
will  then  change  from  compressive  to  tensile. 

Another  example  of  positive  and  negative'  temperature  curvature, 
important  to  the  wind  tunnel  facility,  is  shown  in  Figure  19 .  This 
figure  shows  the  radial  temperature  distribution  in  the  matrix  and 
insulation  of  a  storage  heater.  Note  that  this  curve  has  a  negative 
second  derivative  within  the  matrix,  but  a  positive  second  derivative 
within  the  insulation. 

One  result  of  the  analysis  was  to  determine  the  effect  of 
positive  and  negative  temperature  curvature  on  the  stress  field  in  un¬ 
restrained  blocks.  As  will  be  seen,  positive  curvature  causes  compres¬ 
sive  stresses  at  the  surface,  and  negative  curvature  causes  tensile 
stresses  at  the  surface. 

4.2.2  Block  Shane  and  Size 

Shape  and  size  of  an  elastic  body  can  significantly  affect  the 
magnitude  of  thermal  stress  it  experiences.  In  the  gross  sense,  in¬ 
creases  in  body  dimensions  will  produce  larger  deviations  (from  linear) 
of  the  temperature  profile  and  also  dictate  the  development  of  larger 
internal  strain.  Both  of  these  factors  cause  an  increase  in  stress.  The 
shape  of  a  body  is  important.  If  one  of  the  dimensions  are  small  com¬ 
pared  to  the  other  two,  the  stress  field  will  tend  to  be  two-dimensional. 
A  change  in  shape  can  significantly  affect  the  stress  field,  depending 
on  the  orientation  of  the  change  with  respect  to  the  direction  in  which 
the  temperature  varies . 

For  a  given  one-dimensional  temperature  distribution  and 
a  given  plan  form  of  brick,  the  computer  analyses  investigated  three 
bricks  1"  x  1"  in  planform,  and  1",  1/2",  and  1/4"  thickness.  Con¬ 
sideration  of  the  thickness  variable  is  important  because  of  the  follow¬ 
ing  reasons.  The  best  design  will  result  from  a  compromise  of  all  three 
considerations . 
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1 .  The  temperature  at  the  back-side  (cold  side)  must  be 
low  erlfru&h  to  be  compatible  with  the  adjacent  material, 
and  also  compatible  with  the  type  of  bonding  or  fasten¬ 
ing  to  be  used. 

2.  It  is  desirable  that  strains  on  the  backside  surface  are 

at  a  low  level  to  prevent  bond  failure  between  the  ceramic 
blocks  and  the  backup  material. 

3.  The  thickness  variable  will  have  an  effect  on  the  maximum 
thermal  stress  magnitudes.  The  thickness  should  be  small 
enough  to  reduce  maximum  tensile  and  shear  stresses  (or, 
more  properly,  strains)  to  a  tolerable  level. 

For  a  given  one  dimensional  temperature  distribution  and  a  given 
thickness  of  block,  the  magnitude  and  distribution  of  thermal  stresses 
will  vary  with  a  change  in  the  planform  size  and  shape  of  a  block.  The 
computer  analysis  investigated  a  block  1"  thick  and  2"  x  2"  in  planform 
(hot  face) ,  then  subdivided  the  block  as  shown  in  Figure  20  •  These 
results  provide  information  needed  to  establish  a  maximum  planform 
size  and  shape. 

4.3  COMPUTER  PROGRAM 

The  program  developed  obtains  thermal  stresses  created  within 
a  rectangular  block  by  the  variation  of  temperature  and  material  properties 
in  one  direction .  The  governing  equations  are  formulated  by  means  of 
the  finite -element  method  and  then  solved  by  the  Gauss-Seidel  iterative 
method.  A  description  of  the  theory  and  program  TPA  appears  in  Appendix 
I.  A  revised  program  TPA  II  is  described  in  Appendix  II. 

The  method  essentially  divides  the  block  into  a  number  of  elements 
as  shown  in  Figure  21.  Symmetry  allows  analysis  of  only  one  quadrant 
of  the  block,  permitting  greater  accuracy  for  a  given  number  of  elements. 

The  original  program,  termed  TPA  was  in  FORTRAN  IV  for  the 
IBM  7094  computer.  The  revised  program,  TPA  II,  is  also  coded  in 
FORTRAN  IV  but  for  use  on  the  CDC  6600  computer.  TPA  II  permits  a 
much  greater  subdivision  of  the  blocks  and  correspondlingly  greater 
accuracy.  The  CDC  6600  machine  has  a  larger  core  storage,  permitting 
use  of  more  nodes  plus  more  rapid  execution.  The  program  conversion 
was  performed  by  IITRI. 

4.3.1  Assumptions 

At  least  three  factors  affect  the  accuracy  of  the  thermal  stress 
calculations » 
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1 .  The  number  of  subdivisions  of  the  rectangular  block 
used  in  calculating  the  thermal  stresses  affects  the 
accuracy  of  the  calculations .  In  general,  the  capa¬ 
bilities  of  program  TPA  II  are  sufficient  to  reduce  errors 
to  acceptable  limits  • 

2.  The  analysis  is  an  elastic  analysis,  which  assumes  that 
stress  is  proportional  to  strain.  The  effects  of  creep  and 
plastic  flow  are  not  considered.  For  the  types  of  ceramics 

»  considered  in  the  report,  plastic  flow  occurs  only  at  high 
temperatures  (>  2000 °F).  Also,  because  of  the  steep 
temperature  profile,  only  a  small  layer  of  the  block  is  at 
temperatures  greater  than  2000° F.  Thus  the  analysis  is 
a  reasonable  and  generally  conservative  approximation 
to  the  true  case. 

3.  Material  property  data  for  temperatures  in  excess  of  2000°F 
are  not  readily  available .  The  analysis  uses  estimated 

and  extrapolated  values  for  high  temperatures.  The  property 
data  most  in  doubt  are  the  modulus  of  elasticity  and 
Poisson's  Ratio.  Thus  these  "estimated"  values  also 
contribute  to  the  uncertainty  of  final  calculated  results . 

All  computer  runs  were  made  for  calcia  stabilized  zlrconia 
with  a  porosity  of  0.05.  Coefficient  of  thermal  expansion 
and  modulus  of  elasticity  data  are  shown  in  Figures  22  and 
23.  Due  to  lack  of  published  data  of  Poisson's  ratio  versus 
temperature,  a  constant  value  of  0.29  was  used.  All  of 
the  above  property  values  were  based  on  Reference  2 . 

4.3.2  Problem  Specifications  and  Input  -  Output  Data 

Table  III  gives  a  summary  of  the  runs  for  the  block  thermal  stress 
analyses.  Succeeding  graphs  and  tables  show  material  property  data, 
temperature  distributions,  and  stress  data  from  the  computer  runs.  A 
total  of  16  computer  runs  were  made.  Runs  1  through  8  were  completed 
by  IITRI  on  an  IBM  7094  computer  using  TPA.  Runs  9  through  16  were 
run  on  a  CDC  6600  computer  by  CDC  in  Minneapolis  using  TPA  II. 
Description  of  input  data  is  given  Tables  IV,  V  and  VI  and  in  Appendixes 
I  and  II. 


The  initial  output  is  an  echo  of  all  input  data  •  Following  this , 
the  applied  nodal  forces  and  starting  values  for  the  displacements  are 
printed  out,  and  then  whatever  intermediate  output  the  user  specifies. 

The  final  output  consists  of  the  displacements  and  stresses. 
The  stresses  are  computed  and  printed  at  the  centroids  of  all  elements. 
Program  TPA  also  gives  stresses  at  the  centers  of  all  outside  element 
faces.  TPA  II  gives  only  stresses  at  the  centroids  of  all  elements. 

IITRI  later  advised  that  the  stresses  at  the  outside  faces  should  be  dis¬ 
regarded  in  TPA.  Better  values  are  obtainable  by  extrapolating  from  the 
interior  to  the  surface . 
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The  output  for  each  given  point  includes: 

1.  The  three  orthogonal  normal  stresses. 

2.  The  three  orthogonal  shear  stresses. 

3.  The  three  principal  stresses  and  their  direction  cosines. 

4.  The  "effective"  stress,  defined  as 


s2  -  bxx  -  d)2  +  byy  -  d)2  +  bzz  ~  11)2 

2  2  2 
+  2a  +  2a  +  2a  . 
xz  yz  xy 


d  =  T  bxx+<Iyy+azz) 


The  interpretation  of  the  results  of  the  analysis  requires  some 
care.  The  stresses  in  an  element  vary  linearly  in  the  direction  of  the 
coordinate  axes.  Since  the  stresses  satisfy  equilibrium  only  on  the 
average  it  is  probably  most  meaningful  to  consider  stresses  only  at 
the  center  of  the  elements.  Estimates  of  stresses  at  outside  faces  can 
best  be  made  by  extrapolating  the  results  obtained  at  the  centers  of 
elements  as  one  proceeds  from  the  interior  to  the  surface. 

4.4  CALCULATION  RESULTS 

This  section  contains  results  of  the  computer  analysis .  It  showns 
effects  of  analytically  varying  the  following  parameters . 

1  •  Effect  of  curvature  of  temperature  profile  on  the 
internal  stress  field. 

2.  Variation  of  stresses  with  block  thickness. 

3.  Variation  of  stresses  with  plan  size  and  shape  of  block. 

4.  Variation  of  stresses  with  time  from  start  of  run. 

5.  Principal  stresses  and  maximum  shear  stresses. 

6.  Variation'of  stresses  with  material  properties. 
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It  may  be  noted  here  that  the  calculated  stress  distributions 
were  examined  to  assure  that  they  satisfied  static  equilibrium.  That 
is, forces  acting  on  any  plane,  and  moments  about  any  line  in  the  plane, 
must  be  zero.  That  the  results  satisfy  equilibrium  was  confirmed  by  hand 
calculations  of  a  number  of  cases . 

4.4.1  Effect  of  Curvature  of  Temperature  Profile 

Calculations  were  made  to  establish  the  qualitative  effect  of 
temperature  profile  curvature  (i.e.,  positive  and  negative  second 
derivatives)  on  stress  distribution.  Three  temperature  curves  were 
used  (Figure  24):  one  with  positive  curvature,  which  is  the  40 -second 
curve  from  Figure  18,  a  reversal  of  this  curve,  giving  negative  curvature, 
and  a  linear  distribution . 

The  stress  distributions  are  shown  in  Figure  25.  Consider  first 
the  results  for  the  positive  curvature  temperature  profile  (at  the  top  of 
the  page) .  Compressive  stresses  exist  at  the  surfaces  of  the  block, 
and  tensile  stresses  exist  near  the  center  of  the  block. 

Examination  of  the  stress  distributions  presented  in  later  figures 
indicates  that  this  is  generally  true.  In  some  cases  it  appears  the  stresses 
approach  zero  at  the  side  of  the  block,  and  may  even  become  slightly 
tensile.  But,  in  general,  positive  curvature  of  the  temperature  distribution 
produces  compressive  stresses  on  the  surface  and  tensile  stresses  in 
the  interior  of  a  block . 

The  linear  temperature  distribution  should  result  in  zero  stresses 
(s ee, for  example.  Reference  3) .  The  calculated  stresses  are  not  zero 
because  of  the  approximations  in  the  compution  method,  especially  the 
finiteness  of  the  mesh. 

The  stress  distribution  for  the  negative  curvature  temperature 
profile  is  qualitatively  the  reverse  of  that  for  the  positive  curvature 
profile.  Tensile  stresses  exist  on  and  near  the  surfaces  of  the  block, 
and  compressive  stresses  exist  near  the  center  of  the  block. 

We  note  also  that  the  maximum  stress  levels  occur  in  the  portion 
of  the  block  where  the  second  derivative  of  the  temperature  distribution 
is  the  highest.  Comparison  of  these  and  other  results  has  also  shown 
that  the  maximum  stress  is  roughly  proportional  to  the  departure  of  the 
temperature  curve  from  linearity. 
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In  summary,  the  following  conclusions  are  made: 


Shape  of  Curve 
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1.  Stress  roughly  proportional  to  AT. 

2 .  Maximum  stress  occurs  in  region  of  maximum  second 
derivative  of  temperature  curve. 

4.4.2  Effect  of  Block  Size  and  Shape 

A  general  procedure  for  sizing  a  block  in  the  throat  insulation 
layer  for  acceptable  thermal  stresses  would  involve: 

1.  Select  the  minimum  thickness  of  block  consistent 
with  allowable  backside  temperature. 

2.  Using  the  above  thickness,  select  a  reasonable  planform 
size  and  shape  that  will  result  in  acceptable  thermal 
stresses. 

Thus  both  thickness  and  planform  are  important.  (Planform  is  the  surface 
exposed  to  the  heat  transfer.) 

The  effect  of  thickness  was  examined  by  calculating  stresses  in 
1  x  1-inch  planform  blocks  having  thicknesses  of  1/4,  1/2  and  1-inch. 
These  results  depend  strongly  on  the  temperature  distribution.  In  this 
case  the  10-second  curve  of  Figure  18_was  used,  with  temperature  de¬ 
pendent  properties.  For  the  1/4  and  1/2 -inch  thicknesses,  only  the 
corresponding  portion  of  the  temperature  curve,  measured  from  the  hot 
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face,  was  used.  The  results  are  shown  in  Figure  26.  As  expected, 
the  stresses  decrease  as  the  thickness  is  decreased.  The  maximum 
stresses  are  plotted  in  the  right  hand  portion  of  Figure  27  . 

The  effect  of  planform  size  and  shape  was  calculated  for 
various  blocks  that  were  all  one  inch  thick,  and  therefore  all  had  the 
same  temperature  distribution.  Thus  the  effect  of  temperature  distri¬ 
bution  was  eliminated.  The  10 -second  temperature  curve  of  Figure  18 
was  used,  with  temperature  dependent  properties. 

The  effect  of  planform  size  is  shown  in  Figure  28.  The  effect 
of  planform  shape  is  shown  in  Figure  29.  The  maximum  stresses  in  both 
cases  are  plotted  on  the  left  hand  side  of  Figure  27. 

It  is  clear  from  these  results  that  a  reduction  of  only  one  plan- 
form  dimension  (Figure  29)  is  not  very  effective  in  reducing  maximum 
stresses.  Both  dimensions  must  be  reduced  (Figure  28).  The  shape 
with  the  lowest  stresses  is  rod-like,  with  the  temperature  varying 
along  the  length  of  the  rod.  (The  1/4  x  l/4  x  1-inch  block  in  Figure 
28  is  equivalent  to  a  rod  suddenly  heated  at  one  end.) 

4.4.3  Variation  of  Stresses  with  Time  During  Run 

Figure  30  illustrates  how  the  magnitude  and  distrlbtuion  of 
thermal  stresses  changed  during  a  wind-tunnel  run.  Initially,  a  thin 
section  of  the  hot  surface  has  very  high  compressive  stresses .  As 
the  run  progresses,  the  thermal  stress  field  grows  outward  away  from 
the  hot  surface.  Maximum  tensile  stresses  Increase  and  maximum 
stresses  (tensile  and  compressive)  shift  away  from  the  hot  surface. 

At  a  later  period  in  time,  both  tensile  and  compressive  stresses  seem 
to  decrease.  This  is  a  result  of  the  temperature  curve  becoming  more 
linear.  The  lower  stresses  are,  however,  offset  by  reduced  strengths 
of  the  material  at  higher  temperatures .  The  curves  illustrate  the  complex 
nature  of  the  thermal  stress  history.  It  is  difficult  to  determine  which 
portion  of  the  run  cycle  represents  a  more  severe  condition  for  the  blocks. 

The  stress  levels  are  much  higher  than  the  strengths  of  ceramic 
materials .  In  an  actual  run  the  material  would  crack  and  thereby  would 
relieve  stresses.  The  relation  between  the  calculated  stress  field  and 
spalling  is  as  yet  unknown.  The  early  high  compressive  stresses  also 
indicate  high  shear  stresses  at  the  surface,  which  would  contribute  to 
spalling.  The  intermediate  result  shows  high  tensile  stresses  in  the 
block  interior.  And' finally,  during  the  cool  down  cycle,  tensile  stresses 
will  occur  on  the  outside  surfaces  of  the  block.  These  could  also  con¬ 
tribute  to  cracking  and  spalling. 

4.4.4  Principal  Stresses  and  Maximum  Shear  Stresses 

Figures  31  and  32  show  the  principal  stress  and  maximum  shear 
stresses  for  a  1  x  1  x  1-inch  block  and  the  10 -second  temperature  curve. 
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Information  is  given  for  X  =  0 .07  and  X  =  0 .3,  which  are  at  or  near 
the  planes  of  maximum  compressive  and  tensile  stresses. 

Comparison  of  these  data  with  the  1  x  1  x  1-inch  block  of 
Figure  26  shows  that  the  principal  stresses  are  almost  identical  in 
magnitude  with  the  stresses  plotted  parallel  and  perpendicular  to 
the  sides  of  the  block.  This  means  that  for  the  cases  studied,  the 
principal  stress  vectors  tend  to  have  directions  nearly  the  same  as 
the  three  principal  coordinate  axes . 

The  maximum  shear  stress  is  equal  to 

the  maximum  principal  stress  -  the  minimum  principal  stress 

2 

and  acts  in  the  plane  which  bisects  the  angle  between  the  maximum 
and  minimum  principal  stresses.  The  shear  stress  is  important  because, 
in  areas  where  the  block  is  under  high  compressive  thermal  stresses, 
shear  may  be  the  mode  of  failure.  Close  examination  of  the  direction 
cosines  for  the  principal  stresses,  and  the  corresponding  directions  of 
the  maximum  shearing  stresses,  suggest  that  shear  failures  would  occur 
in  a  manner  to  round  off,  or  break  comers  off,  the  rectangular  block. 

This  is  also  intuitive,  and  in  fact  does  often  happen. 

4.4.5  Variation  of  Stresses  with  Material  Properties 

The  purpose  of  this  comparison  was  to  evaluate  the  effect  of 
reduced  modulus  of  elasticity  at  high  temperatures  on  the  overall  thermal 
stress  pattern.  Both  plots  shown  on  Figure  33  used  the  40-second  tempera¬ 
ture  curve  and  were  calculated  using  nonequal  layer  divisions . 

The  stress  fields  show  that  utilization  of  correct  material  proper¬ 
ties  is  important  for  the  hot  surface  region.  For  the  case  in  question, 
stresses  differ  by  a  factor  of  two  in  the  high  compressive  regions .  How¬ 
ever,  note  that  thermal  stresses  in  the  remainder  of  the  block  are  not 
significantly  different. 

Also  note  that  the  results  based  on  temperature  dependent  material 
properties  have  reduced  thermal  stresses  at  the  surface.  This  arises 
from  the  fact  that  modulus  of  elasticity  decreases  more  than  the  total 
strain  increases  near  the  surface. 

4.4.6  Computer  Program  Use 

As  previously  indicated,  two  variations  of  the  program  are  avail¬ 
able  for  stress  calculations:  TPA  and  TPA  II.  TPA  II  is  the  desired  pro¬ 
gram  both  in  terms  of  capability  and  overall  cost. 
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For  blocks  near  cubic  in  shape,  and  using  5  by  5  nodal  points 
with  10  layers,  and  an  error  tolerance  of  1,500  E-04,  the  program 
converges  in  about  50  to  70  iterations.  Increasing  the  number  of  layers 
to  20  about  doubles  the  number  of  iterations  required  (i.e.,  about  100). 
When  the  element  shape  deviates  significantly  from  a  cubic  shape,  the 
number  of  iterations  becomes  greater,  and  may  take  considerably  more 
than  100  to  converge.  Also,  for  a  given  error  tolerance,  the  number  of 
iterations  required  for  convergence  will  increase  as 'the  total  temperature 
difference  across  the  block  increases  . 

Three  runs  were  made  with  linear  distributions  to  evaluate  the 
calculation  accuracy.  Table  VII  shows  the  results  compared  with  results 
from  the  10-second  temperature  curve. 

This  table  shows  that  a  5  by  5  node  and  10  layer  calculation  may 
give  errors  of  about  5  percent  at  the  center  of  a  block .  At  the  comer 
of  the  block  the  percentage  error  is  much  larger,  about  30  percent,  but 
the  absolute  magnitude  of  the  stress  error  is  larger  by  only  a  factor  of 
two .  Doubling  the  number  of  layers  reduces  the  errors  by  a  factor  of 
two. 


Also  the  table  indicates  that  the  absolute  accuracy  of  the  calcu¬ 
lation  depends  on  the  slope  of  the  temperature  curve.  If  the  temperature 
gradient  is  reduced  by  a  factor  of  two,  the  magnitude  of  error  is  reduced 
by  the  same  factor. 

The  approximate  cost  of  running  a  TPA  II  problem  varies  from  about 
$35  for  a  50  iteration  calculation  to  about  $100  for  a  150  iteration  calcu¬ 
lation. 


4.5  CONCLUSIONS  FROM  STRESS  ANALYSIS 

This  portion  of  the  program  involved  the  preparation  and  utilization 
of  a  computer  program.  The  program  calculates  the  three-dimensional 
thermal  stress  distribution  in  an  unrestrained  rectangular  block  having 
a  one -dimensional  temperature  distribution  and  having  temperature  depen¬ 
dent  properties . 

In  application  to  the  uncooled  throat  design  problem,  the  block 
is  considered  to  be  one  element  of  the  throat  insulation  liner.  The  blocks 
would  be  arranged  in  a  mosaic  pattern  to  form  a  continuous  insulation 
layer. 


A  number  of  calculations  were  made  to  examine  the  influence  of 
the  shape  of  the  temperature  distribution  and  the  size  and  shape  of  the 
block  on  thermal  stress.  The  results  of  these  calculations  are  as  follows: 
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1 .  Temperature  distributions  having  negative  second  derivatives 
produce  tensile  stresses  on  the  surface  of  the  block,  and 
compressive  stresses  near  the  center  of  the  block. 

2 .  Temperature  distributions  having  positive  second  derivatives 
produce  the  reverse  effect,  compressive  stresses  on  the  surface 
and  tensile  stresses  inside  of  the  block. 

3.  Thermal  stresses  are  roughly  proportional  to  the  deviation  of  the 
temperature  distribution  from  linearity  (a  linear  distribution 
produces  no  stresses  if  the  properties  are  uniform). 

4.  Thermal  stresses  can  be  reduced  significantly  by  changes  in 
size  and/or  shape.  Reducing  overall  size  is  a  common  approach 
to  reduction  of  stresses  •  But  if  only  one  or  two  dimensions 

can  be  reduced,  the  orientation  with  respect  to  the  temperature 
distribution  is  important.  The  lowest  stresses  are  produced  in 
a  block  that  has  a  rod -like  shape,  with  the  temperature  varying 
along  the  length  of  the  rod . 

5.  The  calculated  stresses  were  generally  much  larger  than  the 
strength  of  available  high  temperature  materials  •  This  in 
itself  does  not  imply  failure,  but  only  that  the  parts  will  crack. 
Cracks  tend  to  reduce  stresses  but  also  to  produce  spalling. 

The  relationship  between  the  calculated  stress  distributions 
and  spalling  was  not  studied . 
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SECTION  V 
CONCLUSIONS 


The  basic  purpose  of  the  work  reported  herein  was  to  determine 
the  feasibility  of  an  uncooled  throat  for  a  hypersonic  wind  tunnel  facility. 
Specific  results  from  the  test  program  and  from  the  thermal  stress  analysis 
are  presented  in  Section  3.5  and  4.5.  The  conclusions  stated  in  this 
section  apply  to  the  fullscale  throat. 

Tests  were  not  made  at  conditions  as  severe  as  would  be  en¬ 
countered  by  a  completely  uncooled  throat  in  the  fullscale  nozzle . 
Furthermore,  the  results  cannot  be  readily  extrapolated  to  fullscale 
conditions.  Thus,  the  conclusions  that  follow  necessarily  involve 
engineering  judgement. 

On  the  basis  of  this  work,  and  other  experience  with  ceramics, 
it  is  the  authors 1  opinion  that  a  completely  uncooled  throat  for  a  large 
wind  tunnel  is  not  feasible  using  currently  available  commercial  materials  • 
However,  the  results  do  indicate  that  materials  are  available  for  a  ceramic 
lined  throat  that  could  operate  at  the  maximum  tunnel  flow  conditions , 
but  with  the  thermal  shock  conditions  reduced  by  the  use  of  film  cooling 
and  by  preheating  through  flow  of  air  during  heater  pressurization. 

Film  cooling  and  preheating  both  reduce  the  rate  of  heat  transfer 
and,  consequently,  the  severity  of  the  thermal  shock.  Film  cooling  also 
reduces  the  maximum  temperature  reached  by  the  ceramic  insulation. 

The  current  facility  design  concept  provides  for  film  cooling  of  the  hot 
valve  seat  with  a  flow  rate  equal  to  9  percent  of  the  mainstream  flow. 

This  film  cooling  (in  combination  with  preheating)  would  probably  be 
adequate  for  a  ceramic  lined  throat  section.  In  comparison,  the  film 
cooling  required  for  the  backside  cooled  throat  is  25  percent  of  the  main¬ 
stream  flow. 

The  tests  indicate  that  several  zirconia  and  zirconium  diborlde 
compositions  are  available  that  could  withstand  the  flow  conditions  with 
film  cooling  and  preheating.  Application  of  these  materials  to  a  fullscale 
throat  poses  design  problems  that  have  not  been  solved.  The  basic 
concept  of  a  liner  composed  of  many  small  blocks  (nominally  one  inch 
cubes)  involves  attachment  of  the  blocks  to  a  backup  structure.  This 
would  be  the  most  difficult  design  problem  in  applying  the  zirconia 
materials . 

The  tests  indicate  that  several  zirconia  and  zirconium  diborlde 
compositions  are  available  that  could  withstand  the  flow  conditions 
with  film  cooling  and  preheating.  Application  of  these  materials  to  a 
fullscale  throat  poses  design  problems  that  have  not  been  solved.  The 
basic  concept  of  a  liner  composed  of  many  small  blocks  (nominally 
one  inch  cubes)  involves  attachment  of  the  blocks  to  a  backup  structure. 
This  would  be  the  most  difficult  design  problem  in  applying  the  zirconia 
materials . 
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The  zirconium  diboride  materials  have  sufficient  oxidation 
resistance,  if  used  with  film  cooling,  because  of  the  effect  of  film 
cooling  in  reducing  the  temperature  level.  However,  these  materials 
have  high  thermal  conductivity,  equal  to  or  greater  than  carbon  steel  • 

Since  they  are  not  insulators,  their  use  would  require  a  different  design 
concept.  It  may  be  possible  to  use  the  diborlde  materials  in  combination 
with  a  zirconia  insulation  liner. 

A  computer  program  was  prepared  and  used  for  calculation  of 
thermal  stresses  in  blocks  that  would  form  the  insulation  liner.  These 
blocks  are  characterized  by  a  one-dimensional  temperature  distribution 
(normal  to  nozzle  surface)  and  by  having  essentially  no  mechanical 
restraint. 

Calculation  results  provided  a  description  of  the  stress  distribution 
and  the  influence  of  size  and  shape  of  the  blocks  •  Stresses  can  be  re¬ 
duced  most  effectively  by  reducing  all  dimensions  of  the  surface  exposed 
to  the  hot  flow.  That  is,  cumulative  thermal  expansion  should  be  as 
small  as  possible .  Minimum  stresses  are  produced  in  rod  shaped  pieces, 
with  the  rod  axis  parallel  to  the  heat  flow . 
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TABLE  1 

SUMMARY  OF  TEST  SPECIMEN  MATERIALS 
TEST  PHASE  I 

Test 

Specimen  Description 


1  and  2 

Tungsten- Rhenium  wire  reinforced,  lsostatically  pressed 
partially  stabilized  MgO-CaO  zirconia. 

Density:  350  lbs/ft3 .  Supplier:  Thompson  Ramo  Woolridge,  Inc. 

3 

60D  castable,  partially  stabilized  calcla -zirconia. 

Maximum  grain  size  60  mesh. 

Density:  260  lbs/ft3 .  Supplier:  Zircoa 

4 

60 D  castable,  fully  stabilized  yttria-zirconia . 

Maximum  grain  size  60  mesh. 

Density:  260  lbs /ft  3 .  Supplier:  Zircoa 

5  and  6 

Fine  grain,  dense,  pressed  partially  stabilized  magnesla- 
calcia -zirconia,  Zircoa  composition  1027. 

Density:  340  lbs/ft’.  Supplier:  Zircoa 

7  and  8 

8D  castable,  partially  stabilized  calcia-zriconia . 

Maximum  grain  size  8  mesh 

Density:  260  lbs/ft3.  Supplier:  Zircoa 

9  and  10 

8D  castable,  partially  stabilized  magnesia-calcla -zirconia, 

Zircoa  composition  1027 .  Maximum  grain  size  8  mesh. 

Density:  230  lbs/ft3.  Supplier:  Zircoa 

11 

60 D' castable,  partially  stabilized  magnesia-calcia -zirconia, 
Zircoa  composition  1027.  Maximum  grain  size  60  mesh. 

Density:  220  lbs/ft3.  Supplier;  Zircoa 

12  and  13 

Pressed  partially  stabilized  magnesia  zirconia,  Zircoa 
composition  1721. 

Density:  260  lbs/ft3.  Supplier:  Zircoa 

TEST  PHASE  II 

14  and  15 

Pressed  10.4%  yttria-zirconia,  Coors  composition  ZP-997  grogged 
Density:  296  lb/ft3.  Supplier:  Coors  Porcelain  Co. 

16  and  17 

Same  as  14  and  15 

18  and  19 

ManLabs  Diboride  Material  V,  zirconium  dlboride  with 

20  volume  percent  silicon  carbide. 

Density:  340  lb/ft3.  Supplier:  ManLabs,  Inc. 

20  and  21 

Same  as  18  and  19 

22  and  23 

ManLabs  Dlboride  Material  VIII,  zirconium  diboride  with  14 
volume  percent  silicon  carbide  and  30  volume  percent  carbon. 
Density:  330  lb/ft3.  Supplier:  ManLabs,  Inc. 

24  and  25 

Same  as  22  and  23 

32 


TABLE  II. 

SUMMARY  OF  UNCOOLED  THROAT  RUNS 
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PHASE  I 


P 

o 

T 

0 

Press. 

Time 

Run 

Time 

Depress . 
Time 

Specimens 

Run 

psla 

°R 

Sec 

sec 

sec 

Left 

Right 

1 

345 

1800 

30 

0 

130 

1 

2 

2 

400 

1700 

60 

19 

140 

1 

2 

3 

300 

2350 

23 

39 

200 

1 

2 

4 

415 

2200 

28 

26 

280 

3 

4 

5 

300 

2050 

29 

3 

43 

5 

6 

6 

415 

1950 

36 

35 

130 

5 

6 

7 

405 

2050 

24 

34 

152 

5 

6 

8 

415 

2200 

35 

32 

140 

7 

8 

9 

415 

2200 

54 

32 

131 

9 

10 

10 

415 

2800 

40 

31 

140 

9 

10 

11 

415 

3100 

45 

36 

145 

5 

11 

12 

415 

2750 

31 

34 

140 

12 

13 

13 

415 

3300 

37 

34 

220 

12 

13 

14 

415 

3000 

29 

39 

205 

9 

10 

15 

415 

3300 

33 

3 

94 

8 

11 

PHASE  II 


p 

T 

Press . 

Run 

Depress . 

r 

O 

l 

O 

Time 

Time 

Time 

Specimens 

Run 

psia 

°  R 

Sec 

sec 

sec 

Left 

Right 

1 

445 

3035 

42 

31 

148 

14 

15 

2 

325 

3180 

40 

1 

96 

18 

19 

3 

'400 

3160 

52 

32 

141 

18 

19 

4 

400 

3180 

37 

21 

98 

22 

23 

5 

400  . 

3140 

48 

37 

138 

16 

17 

6 

600 

3160 

52 

32 

141 

22 

23 

7 

800 

3260 

62 

31 

147 

22 

23 

8 

700 

3350 

48 

2 

151 

22 

23 

9 

600 

3550 

41 

8 

139 

18 

19 

33 


TABLE  m. 

SPECIFICATIONS  FOR  THERMAL  STRESS  COMPUTER  RUNS 


Run 

Time  (sec) 

Dimensions 

Variables 

Studied 

Comments 

a 

b 

Thick¬ 

ness 

1 

0.1 

1.0 

1.0 

1.0 

T 

2 

10.0 

1.0 

1.0 

1.0 

T,  t,  s 

3 

40.0 

1.0 

1.0 

1.0 

T,  m 

4 

10.0 

2.0 

2.0 

1.0 

s 

5 

10.0 

0.5 

0.5 

1.0 

s 

6 

10.0 

1.0 

1.0 

0.25 

t 

7 

10.0 

1.0 

0.5 

1.0 

s 

8 

40.0 

1.0 

1.0 

1.0 

m 

Same  thickness  division  as  Run  3 

9 

40.0 

1.0 

1.0 

1.0 

c 

Same 't'  dlv.  as  Runs  10  and  11 

10 

See  Comments 

1.0 

1.0 

1.0 

c 

Linear  temp,  curve,  10  layers 

11 

See  Comments 

1.0 

1.0 

1.0 

c 

Inverted  Temp.  Curve 

12 

10.0 

1.0 

1.0 

0.5 

t 

13 

10.0 

0.25 

0.25 

1.0 

s 

14 

10.0 

1.0 

0.25 

1.0 

3 

15 

See  Comments 

1.0 

1.0 

1.0 

d 

Linear  temp,  curve,  20  layers 

16 

See  Comments 

1.0 

1.0 

1.0 

temp  .gradient 

Linear  temp,  curve,  10  layers 

Legend 

T  -  Time  s  =  Plan  form  size  or  shape  m  =  Material  properties 

t  =  Thickness  c  =  Shape  of  temp  .curve  d  =  Subdivisions  of  block 
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TABLE  IV. 

INPUT  DATA  FOR  COMPUTER  RUNS  9,  10  AND  11 


Constant  Properties  Evaluated  at  Room  Temperature 
E  =  26.5xl06  a  -  4.55  x  10~6 

Poisson's  Ratio  =  0.29 


RUN  9 

RUN  10 

RUN  11 

layer 

Thickness 

H 

Temp 

T 

Layer 

Thickness 

H 

Temp 

T 

Layer 

Thickness 

H 

Temp 

T 

1 

0.1 

3250 

1 

0.1 

3730 

1 

0.1 

3870 

2 

0.1 

2130 

2 

0.1 

3345 

2 

0.1 

3780 

3 

0.1 

1350 

3 

0.1 

2960 

3 

0.1 

3660 

4 

0.1 

800 

4 

0.1 

2575 

4 

0.1 

3540 

5 

0.1 

460 

5 

0.1 

2190 

5 

0.1 

3370 

6 

0.1 

260 

6 

0.1 

1805 

6 

0.1 

3200 

7 

0.1 

150 

7 

0.1 

1420 

7 

0.1 

2980 

8 

0.1 

110 

.  8 

•  0.1 

1035 

8 

0.1 

2720 

9 

0.1 

80 

9 

0.1 

650 

9 

0.1 

2350 

10 

0.1 

70 

10 

0.1 

265  ’ 

10 

0.1 

1410 
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TABLE  V. 

INPUT  DATA  FOR  COMPUTER  RUN  12 


Young 1  s  Thermal 

Layer  Thickness  Temperature  Modulus,  Coefficient 


No. 

H 

T 

r 

E 

a 

1 

.020 

/ 

3480 

1.8  x  106 

7.2  x  10“6 

2 

.020 

.  2880 

4.0  x  106 

7.0  x  10“6 

3 

.020 

2440 

10.0  x  106 

6.8  x  10“6 

4 

.020 

2070 

19.3  x  106 

6.6  x  10"6 

5 

.020 

1760 

22.0  x  106 

6.4  x  10“6 

6 

.050 

1260 

24.5  x  106 

6.0  x  10“6 

7 

.050 

740 

25.8  x  106 

5.5  x  10"6 

8 

.100 

300 

26.5  x  106 

5.0  x  10-6 

9 

.100 

150 

26.5  x  106 

4.8  x  10“6 

10 

.100 

80 

26.5  x  106 

4.6  x  10“6 

Time  =  10  sec. 
Thickness  =  0.50  in. 
Plenform  Size  =  1"  x  1" 
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TABLE  VI. 

INPUT  DATA  FOR  COMPUTER  RUNS  15  AND  16 
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Use  Room  Temperature  Properties  for  All  Layers 
E  =  26.5  x  106 
a  =  4.5  x  10“6 

Poisson's  Ratio  =  0.29 


RUN  15 

RUN  16 

Layer 

Thickness 

H 

Temp 

T 

Layer 

Thickness 

H 

Temp 

T 

1 

0.05 

3020 

1 

0.1 

1993 

2 

3627 

2 

1801 

3 

3434 

3 

1609 

4 

3241 

4 

1417 

5 

3048 

5 

1225 

6 

2855 

6 

1033 

7 

2662 

7 

841 

8 

2469 

8 

649 

9 

2276 

9 

457 

10 

2083 

10 

265 

11 

1890 

12 

1697 

13 

1504 

14 

1311 

15 

1118 

16 

925 

P 

732 

18 

539 

19 

346 

20 

153 
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TABLE  VII. 

COMPARISON  OF  CALCULATED  STRESSES  RESULTING  FROM  A  LINEAR 
TEMPERATURE  CURVE  AND  THE  40 -SECOND  TEMPERATURE  CURVE 


Temperature 

Distribution 

Nodes 

and 

Layers 

♦Max. 

Compressive 
Stresses 
Center  of 

Hot  Surface 
(KSI) 

*Max. 

Compressive 
Stresses 
Comer  of 
Hot  Surface 
(KSI) 

40  sec.  curve 

5x5x10  layers 

105 

36 

linear-same  slope 

5  x  5  x  10  layers 

6 

11 

linear-same  slope 

5  x  5  x  20  layers 

3 

7 

linear-half  slope 

5  x  5  x  10  layers 

3 

5 

♦Stresses  parallel  to  hot  surface 
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FIGURE  1 .  FACILITY  SETUP  FOR  UNCOOLED  THROAT  TESTS 
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FIGURE  2 
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FIGURE  3.  TEST  APPARATUS  COMPONENTS 
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FIGURE  4 
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FIGURE  7 
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00 


i  A 


Test  Specimen  No.  1 
and  2  (both  identical) 


TUNGSTEN  REINFORCED  ZIRCONIA 


Test  Specimen  No.  A 


60  D  CASTABLE  YTTRIA  ZIRCONIA 

FIGURE  8.  PRETEST  MATERIAL 


Test  Specimen  No.  3 


60  D  CASTABLE  ZIRCONIA 


FINE  GRAIN  DENSE  1027  ZIRCONIA 


PHOTOGRAPHS 
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CD 


Test  Specimen  No.  11 

60 D  CASTABLE  1027  ZrRCONIA 


Test  Specimen  No.  12 
and  13  (both  identical) 


PRESSED  1721  ZIRCONIA 


8D  CASTABLE  1027  ZIRCONIA 


Test  Specimen  No.  7 
and  8  (both  identical) 


Test  Specimen  No.  9 
and  10  (both  Identical) 


8D  CASTABLE  ZIRCONIA 


FIGURE  9  .  PRETEST  MATERIAL  SPECIMEN  PHOTOGRAPHS 
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Zirconium  Diboride  Specimens  Manufactured  by  ManLabs 


Yt t r i a-Z i rcon i a  Specimens  Manufactured  by  Coors 
FIGURE  10  .  PRETEST  MATERIAL  SPECIMENS,  PHASE  II 
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TEST  SPECIMENS  1  AND  2 
AFTER  RUN  3 


TEST  SPECIMENS  3  AND  4 
AFTER  RUN  4 


TEST  SPECIMENS  5  AND  6 
AFTER  RUN  7 


TEST  SPECIMENS  7  AND  8 
AFTER  .RUN  8 


FIGURE  11.  AFTER  RUN  MATERIAL  SPECIMEN  PHOTOGRAPHS ,  PHASE  I 
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TEST  SPECIMENS  9  AND  10 
AFTER  RUN  10 


TEST  SPECIMENS  5  AND  11 
AFTER  RUN  11 


TEST  SPECIMENS  9  AND  10 
AFTER  RUN  14 


TEST  SPECIMENS  12  AND  13 
AFTER  RUN  13 


FIGURE  12.  AFTER  RUN  MATERIAL  SPECIMEN  PHOTOGRAPHS,  PHASE  I 
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AFTER  RUN  15 
SPECIMENS  8  AND  11 
LOOKING  UPSTREAM 


FIGURE  13.  AFTER  RUN  MATERIAL  SPECIMEN  PHOTOGRAPHS,  PHASE  I 
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Plan  View  of  Test  Specimen  14 


FIGURE  14.  AFTER  RUN  MATERIAL  SPECIMEN  PHOTOGRAPHS,  PHASE  II 
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FLOW 


Test  Specimens  14  and  15  in  Holder 
Looking  Upstream  After*Run  1 


ft 

m  I 

It 

FIGURE  16.  AFTER  RUN  MATERIAL  SPECIMEN  PHOTOGRAPHS,  PHASE  II 
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Test  Specimens  22 
and  23  After  Run  8 


Test  Specimens  18 
and  19  After  Run  9 


FIGURE  17  .  AFTER  RUN  MATERIAL  SPECIMEN  PHOTOGRAPHS,  PHASE  II 
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FIGURE  18.  TRANSIENT  TEMPERATURES  IN  UNCOOLED  THROAT  CERAMIC  ELEMENTS 
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INSULATION 


FIGURE  19.  TEMPERATURE  DISTRIBUTION  IN  THE  BED  MATRIX 

AND  INSULATION  OF  A  CYLINDRICAL  CERAMIC  STORAGE 
HEATER 
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FIGURE  20 .  SUBDIVISIONS  OF  A  2"  x  2"  PLAN  FORM  INTO  SMALLER 
ELEMENTS  TO  STUDY  THE  EFFECT  OF  PLAN-FORM  SIZE 
AND  SHAPE  ON  THERMAL  STRESSES 
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Axis 
of  Block 

XI  Up  To 


ORIGINAL  PROGRAM  (TPA) 


X 


REVISED  PROGRAM  (TPA  II) 


FIGURE  2 1 .  BLOCK  GEOMETRY  AND  SUBDIVISIONS 
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linear  Coefficient  of  Thermal  Expansion  (xlO°)  in. /in.  °F 


i 

a 

FIGURE  22 .  COEFFICIENT  OF  THERMAL  EXPANSION,  a  VS  TEMPERATURE  T 
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FIGURE  23.MODULUS  OF 


Calcia  Stabilized  Zirconia 
Porosity  =  0 .05 


°F 


_ I _ i 

3000 


,  E  VS  TEMPERATURE  T 
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FIGURE  24.  PIECEWISE  CONSTANT  TEMPERATURE  DISTRIBUTIONS 

FOR  STUDYING  EFFECT  OF  KIND  OF  TEMPERATURE  CURVE 
ON  THERMAL  STRESS  DISTRIBUTION  (Computer' Rbns  9  and  10  - 
'CONSTANT  MATERIAL  PROPERTIES) . 
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STRESSES  IN  KSI 
+  TENSION 
-  COMPRESSION 


FIGURE  25.  THERMAL  STRESS  VS.  CURVATURE  OF  TEMPERATURE  PROFILE 
COMPUTER  RUNS  9,  10  AND  11 
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STRESSES  IN  KSI 
+  TENSION 
-  COMPRESSION 


FIGURE  26.  THERMAL  STRESS  VS.  fiLOCK  THICKNESS 
COMPUTER  RUNS  2,  12  AND  6 


Plan  Form  Dimension  "Y"  Height,  Inches 


I nches 


For  10  sec.  Temp.  Curve 


Figure  27 


VARIATION  OP  MAXIMUM  STRESSES  WITH  PLAN  FORM  DIMENSION  &  THICKNESS 
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2X2X1 


ix4xl 


STRESSES  IN  KSI 
+  TENSION 
“  COMPRESSION 


FIGURE  28.  THERMAL  STRESS  VS.  PLAN  FORM  SIZE 
COMPUTER  RUNS  4,  2,  5  AND  13 
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1X£X1 


STRESSES  IN  KSI 
+  TENSION 
-  COMPRESSION 

FIGURE  29.  THERMAL  STRESSES  VS.  PLAN  FORM  SHAPE 
COMPUTER  RUNS  2,  7  AND  14 
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32 


76 


53 


FIGURE  31.  PRINCIPLE  STRESSES  NEAR  HEATED  SURFACE 
COMPUTER  RUN  2 
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FIGURE  32 .  PRINCIPLE  STRESSES  0 .3  -INCHES  FROM  HEATED  SURFACE 
COMPUTER  BUN  2 
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-100 

VARYING  MAT'L 
PROPERTIES 


-96 

CONSTANT  HAT'L 
PROPERTIES 


STRESSES  IN  KSI 
+  TENSION 
-  COMPRESSION 

FIGURE  33.  THERMAL  STRESS  VS.  MATERIAL  PROPERTIES 
COMPUTER  RUNS  3  AND  8 
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FOREWORD 


This  study  o£  the  thermal  stresses  in  rectangular  blocks 
was  performed  by  M.  A.  Salmon  and  T.  Belytschko  of  the 
Structures  Section,  Mechanics  Research  Division  of  IIT  Research 
Institute  (IITRI)  for  the  FluiDyne  Engineering  Corporation. 

Mr.  P.  B.  Hasselquist  coordinated  the  work  for  FluiDyne.  The 

work  was  performed  in  the  period  of  June  to  December  1967. 
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ANALYSIS  OF  THERMAL  STRESSES 
IN  RECTANGULAR  BLOCKS 

ABSTRACT 

An  analysis  of  the  thermal  stresses  in  a  rectangular  block 
by  the  finite  element  method  is  described  and  a  computer  pro¬ 
gram  for  its  implementation  is  presented.  The  program  is 
capable  of  handling  a  temperature  variation  in  one  direction. 
Temperature  dependence  of  the  material  properties  is  taken  into 
account.  Up  to  10  divisions  in  the  direction  of  temperature 
variation  can  be  used;  the  lengths  of  these  divisions  are  speci- 
fied  by  the  user.  Up  to  five  equal  divisions  in  the  other 
two  directions  can  be  used.  Eight  problems  specified  by 
FluiDyne  were  solved  using  the  program. 
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X.  INTRODUCTION 

This  report  describes  a  computer  program  for  the  thermal 
stress  analysis  of  rectangular  blocks  by  the  finite  element 
method.  This  computer  program  has  been  developed  by  IITRI  for 
the  FluiDyne  Corporation  and  applied  to  the  solution  of  eight 
specific  problems.  The  computer  printout  for  these  problems 
has  been  delivered  to  FluiDyne. 

The  method  of  analysis  and  the  computer  program  used  to 
implement  it  are  described  in  the  following  report.  Instruc¬ 
tions  for  the  use  of  the  computer  program  and  the  program  list¬ 
ing  are  given  in  Appendix  A.  Appendix  B  contains  the  data 
specified  by  FluiDyne  for  the  eight  problems  and  the  input 
data  prepared  from  these  specifications  by  IITRI. 

II.  ANALYSIS  AND  COMPUTER  PROGRAM 

An  element  stiffness  matrix  for  a  rectangular  prism  given 
by  Melosh1  is  used  in  the  analysis.  The  prism  geometry  is 
shown  in  Fig.  1.  The  following  expression  Is  used  to  relate 
the  components  of  displacement  to  the  nodal  displacements 


u  ■  4>  X 


(1) 


where 


u  (x,y,z) 


[u  . u  ,u 
X*  y *  zj 


X 


u 


(1) 


.  U. 


X  - “X  * 


(8)  u(l)  ,,(8).  „<1) 


•  .u. 


u. 


*  U, 


(8) 


and 


It  — 

6 

<J)  - 

i 

Hlelosh,  R.  J.,  Structural  Analysis  of  Solids  Proc.  ASCE, 
Vol.  9,  No,  ST4,  August  1963. 
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in  which 

^  -  [-(5+1)  (1-1)  (C+l) ,  (5+1)  (t)+D  (5+D , 
-(5+D  (i+i)  (5-D ,  (5+i)  (i-D  (5-D , 
(5-1)  (rpl)  (5+1)  ,-(5-l)  (l+l)  (5+1) , 
(5-D  (1+1)  (5-D ,-  (5-1)  (i-l)  (5-D  ] 


where 

£  =  x/a,  tj  =  y/b,  C  ■  z/c  . 


The  displacement  field  given  by  Eq.  (1)  satisfies  the  require¬ 
ments  that  displacements  of  contiguous  elements  be  compatible, 
that  rigid  body  displacements  produce  no  strain,  and  that 
states  of  uniform  strain  can  be  represented.  Thus  convergence 
of  the  solutions  to  the  exact  one  as  the  fineness  of  the  mesh 
is  increased  is  assured. 

The  strain-displacement  relations  are 
e  -  D  u  (2) 


where 


£ 


t 


and 


exx* 

eyy»  ezz>  7xy*  7xz* 

7yz  J 

d/dx 

0  0  d/dy 

d/dz 

0 

0 

d/dy  0  d/dx 

0 

d/dz 

0 

0  d/dz  0 

d/dx 

d/dy 

The  stress  strain  law  for  an  isotropic  elastic  material  is 


0  -  C  e 

e 


(3) 
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where  the  elastic  strains  are  expressed  in  terms  of  the  total 
strains  and  the  thermal  strains  by 


e 


+  s. 


and  where 


in  which 


(A) 


c  EO-V) 

*  (l+v)(l-2v) 

Ev 

C2  "  (1+v) (1-2V) 


c  _ _ s _ 

2  (1+v)  * 


The  stresses  derived  from  the  assumed  displacement  field 
of  Eq.  (1)  do  not  satisfy  the  equations  of  equilibrium  in  the 
element.  The  finite  element  method  consists  in  the  selection 
of  a  nodal  force  vector  F  whose  components  correspond  to  the 
elements  of  the  nodal  displacement  vector  X  by  the  application 
of  the  principle  of  virtual  work.  That  is,  the  requirement  is 
enforced  that 


(5) 
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where  X  and  £  are  the  virtual  nodal  displacement  vector  and  the 
corresponding1 strains ,  respectively.  Substitution  of  the  ex¬ 
pressions  for  stress  and  strain  as  functions  of  nodal  displace¬ 
ment  in  Eq.  (5)  gives  the  result. 


F 


DOX-C.E 


(6) 


where  use  has  been  made  of  the  fact  that  the  elements  of  X  may 
be  chosen  in  an  arbitrary  manner. 

Equation  (6)  may  be  written  in  the  form 


K  X  -  F  +  Fp 


(7) 


where 


K 


and 


Dt 


CD® 


V 


JV  d£  c  et 

V 


(8) 


(9) 


are  the  element  stiffness  matrix  and  the  thermal  force  vector, 
respectively. 

The  thermal  strains  are  simply 

ett  -«[l.  °>  0 

where  a  is  the  coefficient  of  expansion  and  T  is  the  tempera¬ 
ture  rise. 

The  formation  of  the  element  stiffness  matrix  in  the  com¬ 
puter  program  is  accomplished  by  the  evaluation  of  Eq.  (7) 

(the  exact  expression  for  the  integrals  involved  are  used) . 

The  equilibrium  equations  for  the  assembly  of  prismatic 
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elements,  into  which  the  block  is  divided,  have  the  same  form 
as  those  for  a  single  element  given  by  Eq.  (7).  The  subdivision 
of  the  block  into  prismatic  elements  is  shown  in  Fig.  2.  The 
program  is  designed  to  handle  problems  in  which  the  temperature 
and  material  constants  (E,  a,  v)  vary  through  the  thickness  of 
the  block  (in  the  x  direction),  but  are  constant  in  any  y-z 

t 

plane.  Owing  to  the  symmetry  of  the  problem  it  is  necessary  to 
consider  only  one  quadrant  of  the  block  as  shown.  Up  to  10 
layers  of  elements  in  the  x  direction  can  be  handled.  The  thick¬ 
ness  of  the  layers  is  specified  by  the  user  as  input.  Up  to 
five  equal  divisions  in  the'x  and  z  directions  can  be  used. 

If  the  maximum  number  of  divisions  is  used  the  block  to  be  an¬ 
alyzed  is  divided  into  250  elements  with  a  total  of  396  nodes. 
The  total  number  of  nodal  displacements  is  1188  of  which  133 
are  known  to  be  zero  (the  displacements  normal  to  the  x-y  and 
x-z  coordinate  planes) .  In  addition  it  is  necessary  to  fix  one 
node  on  the  x  axis  against  displacement  in  the  x  direction. 

In  order  to  accomplish  the  solution  of  this  large  number 
of  equations,  in-core  Gauss  Seidel  iteration  is  used.  A  fur¬ 
ther  reduction  in  storage  requirements  is  obtained  by  taking 
advantage  of  the  symmetry  of  the  stiffness  matrix.  Finally  the 
fact  that  the  elements  in  a  given  layer  have  identical  dimen¬ 
sions  and  material  properties  permits  the  stiffness  coefficients 
to  be  stored  in  a  compact  form.  The  method  used  is  described 
below. 

The  stiffness  matrix  for  a  block  of  n  layers  of  nodes  with 
m  nodes  in  each  layer  can  be  written  in  the  form 


K 


A1  B1 


A2  B2 


SYM. 


An-1  Bn-1 
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Block  Dimensions:  a  x  b  x  c 


Fig. 1-2  SUBDIVISION  OF  QUADRANT  OF  BLOCK  INTO  PRISMATIC  ELEMENTS 
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where  the  A  and  B  are  square  matrices  of  order  3m.  The  elements 
of  the  A  matrices  are  the  forces  on  the  nodes  in  that  layer  due 
to  displacements  of  those  nodes,  while  the  elements  of  the  B 
matrices  are  forces  due  to  displacements  of  nodes  in  the  layer 
above.  Due  to  the  fact  that  the  elements  in  a  layer  are  iden¬ 
tical,  any  node  in  a  layer  can  be  identified  with  one  of  the 
nine  nodes  shown  in  the  sketch. 
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The  A  and  B  matrices  for  this  arrangement  are  of  order  27,  thus 
1458  storage  locations  are  required  for  each  layer. 

The  interpretation  of  the  results  of  the  analysis  requires 
some  care.  The  stresses  in  an  element  vary  linearly  in  the 
direction  of  the  coordinate  axes .  Since  the  stresses  satisfy 
equilibrium  only  on  the  average  it  is  probably  not  meaningful 
to  consider  stresses  anywhere  but  at  the  center  of  the  element. 
Estimates  of  stresses  at  outside  faces  can  best  be  made  by 
extrapolating  the  results  obtained  at  the  centers  of  elements 
as  one  proceeds  from  the  interior  to  the  surface.  The  computer 
program  does  give  stresses  calculated  at  the  centers  of  the  out¬ 
side  faces  of  elements.  However,  examination  of  the  results 
shows  that  these  values  are  not  as  good  as  those  given  by  ex¬ 
trapolation  of  values  at  the  centers  of  elements.  Hindsight 
shows  that  it  was  not  a  good  idea  to  printout  stresses  at  the 
outside  faces  and  these  results  should  be  ignored. 
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APPENDIX  A 
COMPUTER  PROGRAM 

A.  INTRODUCTION 

TPA  is  a  FORTRAN  IV  program  for  obtaining  the  thermal 
stresses  created  within  a  rectangular  prism  by  the  variation  of 
temperature  and  material  properties  in  one  direction.  The 
governing  equations  are  formulated  by  means  of  the  finite  ele¬ 
ment  method  and  then  solved  by  the  Gauss-Seidel  method.  The 
program's  capacity  limits  problems  to  10  layers  of  elements  in 
the  direction  of  the  thermal  gradient  and  up  to  five  elements 
in  each  of  the  other  two  directions. 

B.  INPUT  DATA 

Let  the  direction  of  the  thermal  gradient  be  denoted  by  x 
and  let  the  remaining  two  orthogonal  directions  be  denoted  by 
y  and  z  respectively.  The  x-direction  must  correspond  to  an 

V 

axis  of  the  prism.  Since  the  prism  is  symmetric  with  respect 
to  the  x-y  and  x-z  planes,  only  a  quadrant  of  the  prism  is  con¬ 
sidered.  The  geometry  of  this  quadrant  is  specified  by  giving 

its  y  and  z  dimensions  and  the  number  of  nodes  desired  in  x,  y, 

* 

and  z  directions.  The  nodes  in  the  y  and  z  directions  will  be 
equally  spaced.  The  spacing  between  the  nodes  in  the  x-direction 
which  corresponds  to  the  thicknesses  of  the  element  layers,  are 
specified  by  the  user  in  the  layer  cards. 

For  each  layer  of  elements  the  user  must  provide  a  "layer" 
card  which  gives  its  thickness  and  the  temperature,  Young's 
modulus  E,  Poisson’s  ratio  V,  and  the  thermal  expansion  coeffi- 
c ient  a . 

The  only  other  required  input  is  composed  of  two  control 
numbers  NUREAD  and  NUPRNT  and  the  following  parameters  for  the 
Gauss-Seidel  method:'  the  maximum  number  of  iterations  NIT, 
the  error  printout  cycle  NPERR,  displacement  and  modified  force 
printout  cycle  NPOTPT,  the  tolerance  TOL  and  initial  relaxation 
factor  XFAC.  The  initial  relaxation  factor  should  be  about  1.5. 
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After  the  iterative  method  is  started  if  a  .given  iteration  does 
not  reduce  the  error,  the  relaxation  factor  is  reduced.  The 
Gauss-Seidel  method  terminates  after  the  error  fails  to  decrease 
when  the  number  of  iterations  exceeds  NIT,  or  if  the  error  meets 
the  tolerance  TOL.  During  the  iterations  the  error  and  displace¬ 
ments  are  printed  out  according  to  the  value  of  NPERR  and  NPOTPT, 
respectively.  It  is  recommended  that  the  displacements  not  be 
printed  out  more  than  once  or  twice  during  the  solution,  for 
it  is  very  time-consuming. 

The  control  number  NUREAD  determines  whether  the  starting 
values  of  the  displacements  are  to  be  read  off  a  tape  or  to  be 
computed.  The  second  control  number,  NUPRNT  gives  .the  user  the 
option  of  saving  the  displacements  on  tape  or  interpolating  to 
a  lattice  with  a  larger  number  of  elements.  These  options  allow 
the  user  to  run  the  problem  to  a  certain  point,  save  the  dis¬ 
placements  on  a  tape  and  continue  later  from  the  same  point  by 
using  these  displacements  as  starting  values. 

The  interpolation  routine  inserts  the  displacement  values 
for  a  lattice  with  n-^  elements  in  the  z-direction,  ^  elements 
in  the  y-direction,  and  m  elements  in  the  x-direction  to  obtain 
the  starting  values  for  2n-^  by  2^  by  m  lattice.  When  this 
option  is  used  the  second  problem  is  then  automatically  run 
with  the  interpolated  displacements  as 'starting  values.  A  de¬ 
tailed  description  of  the  input  data  required  is  given  in 
Table  1. 

C.  DESCRIPTION  OF  OUTPUT 

The  initial  output  is  an  echo  of  all  data.  Following  this, 
the  applied  nodal  forces  and  starting  values  for  the  displace¬ 
ments  are  printed  out,  and  then  whatever  intermediate  output  the 
user  specifies. 

The  final  output  consists  of  the  displacements  and  stresses. 
The  stresses  are  computed  and  printed  at. the  centroids  of  all 
elements  and  at  the  centers  of  all  outside  faces. 
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TABLE  1-1 

INPUT  DATA  FORMAT 

Card  1  TITLE  CARD  (12A6) 

Column 

1  ■  72  Any  alphanumeric  information 

Card  2  PARAMETER  CARD  (615,  4F10.0) 

Column 

1  -  5  NX  -  Number  of  nodes  in  x-direction 

6-10  NY  -  Number  of  node6  in  y-airection 

11  -  15  NZ  -  Number  of  nodes  in  z-direction 

16  -  20  NIT  -  Maximum  number  of  iterations 

21  -  25  NPERR  -  Error  printout  cycle 

26  -  30  NPOTPT  -  Displacement  and  modified  force  printout  cycle 
31  -  40  DIMY  -  y-dimension  of  prism  quadrant 

41  -  50  DIMZ  -  z-dimension  of  prism  quadrant 

51  -  60  TOL  -  Error  tolerance  for  Gauss-Seldel 

61  -  70  XFAC  -  Initial  relaxation  factor 

Card  3  LAYER  PROPERTY  CARDS  (15,  5E10.0) 

Column 

1-5  Layer  number 


6-15 

H 

-  Thickness  of  layer 

16  -  25 

T 

-  Temperature  of  layer  / 

26  -  35 

E 

-  Young's  modulus  for  layer 

36  -  45 

PR 

-  Poisson's  ratio  for  layer 

46  -  55 

ALFA 

-  Thermal  expansion  coefficient  for  layer 

Card  4 

CONTROL  CARD 

Column 

1-5 

NUREAD 

-  Determines  how  initial  displacements  are  obtained 

■  0  calculated 

■  1  read  from  card  reader 

■  2  read  from  tape  8 
*  3  read  from  tape  10 

6-10  NUPRNT  -  Gives  the  following  options 

“  0  none  of  the  options  used 

■  1  displacements  are  saved  on  tape  10 

■  2  displacements  are  used  to  interpolate 

to  a  larger  lattice 

-  3  interpolation  and  displacements  from  larger 
problem  saved  on  tape  10 

Card  5  INTERPOLATION  CARD  (OPTIONAL),  (315,  2E10.0) 

Column 

1-5  NIT  -  Maximum  number  of  iterations  for  second  problem 

6-10  NPERR  -  Error  printout  cycle  for  second  problem 

11  u  15  NPOTPT  -  Displacement  printout  cycle  for  second  problem 

16  -  25  TOL  -  Error  tolerance  for  second  problem 

26  -  35  XFAC  -  Initial  relaxation  factor  for  second  problem 
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The  stresses  are  printed  out  in  the  following  manner.  For  each 
point  there  are  five  lines.  The  first  line  gives  the  x,  y  and  z 
coordinates  of  the ‘point.  The  second  line  gives  the  xx,  yy,  zz, 
xy,  xz  and  yz  components  of  the  stresses  and  the  effective' stress 
a  ,  The  effective  stress  a  is  computed  by  the  formula 

CT^  “  (c  *,S)^+(cr  -S)^-t-2(cr  +  CT  +0  ) 

e  v  xx  '  v  yy  *  v  zz  *  v  xy  yz  xz7 

where 

S  *=•  1/3  (cr  +  a  +  c  )  . 
v  xx  yy  zz7 

The  remaining  three  lines  give  the  three  principal  stresses  and 
the  corresponding  principal  directions.  If  the  solutions  to 
additional  problems  are  desired,  the  data  sets  are  repeated. 

Any  number  of  problems  can  be  solved  in  sequence. 
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••••#••»»««•••••»•»••*••**•• . . . .  TPA00020 

THERMAL  STRESS  ANALYSIS  OF  RECTANGULAR  PRISM  TPA00030 

TAPE  8  IS  USED  T-0  INPUT  DISPLACEMENTS  IF 
NUREAO  =  2  OPTION  IS  USED 

TAPE  10  IS  USED  TO  SAVE  DISPLACEMENTS  IF 
NPRINT  =  2  OPTION  IS  USED 


co 

tb 


»•••••••••••• . . . .  .  .  *******JPA00040 

COMMON  A(18, 9,189) 

COMMON / SOL /F ( 21*10,30) , U <21 . 1 0. 30) , I MAX , JMAX , KMAX , SUM, XFAC , TEMP ( 20 
l) » ALFA (20) >H(20) * E I  20 1 »  PR ( 20 ) *  BDI M»  CD  I M 
COMMON/ I DEN /M IDE N( 10)«NIDEN110) .NITER 
DIMENSION  FTEMP16300) ,FTEMPl!6300) .TITLE! 12) 

DOUBLE  PRECISION  F.FTEMP 

EQUIVALENCE (F (1, 1,1) .FTEMP ( l ) ) , CFTEMP1 ( 1 ) .TITLE! 1) ) 


****  READ  AND  PRINT  OATA  **** 

1  READ! 5,100)  T ITLE , I  MAX , JMAX .KMAX ,NI T .NPERR.NPOTPT ,DIMJ. DIMK, TOL , 
1XFAC.DELERR 

IF  (IMAX.EQ.O)  STOP 

100  FORMAT  ! 12A6/6I 5.3F10.0.2F5.0) 

NLAYER-IMAX-1 

READ! 5, 101) !K,H!K) .TEMP! K) ,E ! K ) , PR ! K ) , ALFA l K ) , J= 1 , NL AYER ) 

101  FORMAT! 15, 5FIO.O) 

DIMI=0. 

DO  2  1*1 .NLAYER 

2  DIMI =D IM I £H ( I  ) 

WRITE (6, 200)  TITLE, I  MAX, 01 MI , JMAX.DIMJ, KMAX, DIMK. NIT, NPERR.NPOTPT 


TPAOOHO 

TPA00120 

TPA00130 

TPA00140 

TPA00150 

TPA00160 

TPA00180 

TPA00190 

TPA00200 

TPA00210 

TPA00220 

TPA00230 

TPA00240 

TPA00250 


1 

.  TOL , : 

XFAC.DELERR 

TPA00260 

200  i 

FORMAT! 1H1  12A6/ 

TPA00270 

1 

14H0 

-X-DIRECTION, I4,3X5HN0DES5X7HLENGTH=F8.4/ 

TPA00280 

2 

14H0 

Z-OIRECTION, I4,3X5HN0DES5X7HLENGTH*F8.4/ 

TPA00290 

3 

14H0 

Y-DIRECT I  ON, 1 4 ,3X5HNOOES5X7HLENGTH=F8 .4/ 

TPA00300 

4 

30H0 

MAXIMUM  NO.  OF  ITERATIONS - 13/ 

TPA00310 

5 

30H0 

ERROR  PRINT  CYCLE - 13/ 

TPA00320 

6 

30H0 

OUTPUT  PRINT  CYCLE - 13/ 

TPA00330 

7 

30H0 

ERROR  TOLERANCE - E12.4/ 

TPA00340 

8 

30HO 

OVER  RELAXATION  FACTOR - F7.3/ 

9 

30H0 

DELERR  - F7.3) 

WRITE! 6, 201 ) 

TPA00360 

201  FORMAT  (1H040X18H  LAYER  PRQPERTIES/11H0 
11 1HTEMPERATURE1 3X7HMODULUS5X15HPOI SSON  S 


LAYER  NO. 1 1X9HTHICKNESS9XTPA00370 
RAT I04X 16HEXPANS ION  C0EFFTPA00380 
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n  n  o  n  n  o  ooo  non 


2.  ) 

WRITE (6,203 ) ( K,H(K) , TEMPI K) ,E ( K ) , PR ( K ) , ALFA(K)»K=1,NLAYER) 

203  FORMAT  ( 18* 3X ,5E20. 5) 

300  READ  (5,204)  NUREA0, NUPRNT 

204  FORMAT  ( 1015) 

WRITE  (6,205)  NUREAO, NUPRNT 

205  FORMAT  ( 1HO,5X,8HNUREAO  *,12,5X,8HNUPRNT  *,I2) 

#***  FORM  EQUATIONS  **•* 

XNJ*  JMAX-1 
XNK-  KMAX-1 
BOIM=DIMK/XNK 
CDIM=DIM J/XNJ 
CALL  FSEQ 
NITER«0 

NUREAD-NUREA0S1 

GO  TO  (60,60,80,80) .NUREAD 

••••  MAKE  INITIAL  ESTIMATE  OF  DISPLACEMENTS  **** 

60  CALL  DI$PL(H, TEMP, ALFA, 0 1 M J,DIMK, I  MAX , JMAX ,KMAX ,U,DIMI ) 

GO  TO  91 

80  READ  (8)  (( (U(IfJfK) ,K=1,30J , J=1 , 10 ) , I* 1 , IMAX ) 

91  CONTINUE 

****  ITERATION  ON  NODAL  DISPLACEMENTS  **** 

SF»0. 

DO  250  J-1,6300 
FTEMPK  J)=F  TEMPI  J) 

250  SF  =  SF  E  DABS(FTEMP(J) ) 

SUMP*l.EC35 

DELFAC*  (XFAC-0.8)/10.0 
K0UNT=0 

IF  (NIT.EQ.O)  GO  TO  401 

SET  UP  ROW  AND  COLUMN  IDENTIFIERS 

M IDEN ( 1 ) —  1 
DO  '20  M=2» JMAX 
20  M IDEN( M ) *2 
MIDEN( JMAX)=3 
NIDEN( 1 ) =1 
DO  21  N=2,KMAX 


TPA00390 

TPA00400 

TPA00410 

TPA00420 

TPA00430 

TPA00440 

TPA00450 

TPA00460 

TPA00470 

TPA00480 

TPA00490 

TPA00500 

TPA00510 

TPA00520 

TPA00530 

TPA00540 

TPA00550 

TPA00570 

TPA00580 

TPA00590 

TPA00600 

TPA00650 

TPA00690 

TPA00700 

TPA00710 

TPA00720 

TPA00730 


TPA00750 

TPA00760 

TPA00770 

TPA00780 

TPA00810 

TPA00830 

TPA00840 

TPA00850 

TPA00860 

TPA00870 

TPA00880 

TPA00890 

TPA00900 

TPA00910 


AEDC-TR-70-92 


no  non 


21  N IDEN ( N )  =2 
NIOEN(KMAX) *3 

10  SUM=0. 0 

NI TER=N I TER61 
00  11  L  1  =1 » IMAX 
00  11  Ml=l» JMAX 
DO  11  N 1  =1 »  KM AX 
I3=1£(L1-1)*9 

CALL  SOL VE(A( 1,1,13) , LI, Ml, Nl) 

11  CONTINUE 
IMID=IMAX/2 
UMID=U( IMID,1 ,1) 

DO  1101  Ll=l, IMAX 
DO  1101  M1=1,JMAX 
DO  1101  N 1= 1 , KMAX 

1101  U(L1,M1,N1)=U(L1,M1,N1)-UMID 
DO  601  J=1 » 6300 
601  FTEMP ( J ) =FTEMP1 ( J ) 

****  CYCLE  COUNT  AND  PRINT  CHECK  *•** 

SUM=SUM/SF 

CHECK  WHETHER  RELAXATION  FACTOR  IS  TO  BE  MODIFIED 

IF  (  ABS(SUM) .LT.l.E-10  )  GO  TO  12 
IF  (  (  (SUMP-SUMJ/SUM)  .LT.DELERR  >'  GO  TO  55 
K0UNT*K0UNT£1 
IF  ( KOUNT.GE. 10)  GO  TO  52 
XFAC-XFAC-DELFAC 
GO  TO  55 

52  IF  ( SUM. GT. SUMP)  GO  TO  12 
55  SUMP=SUM 

IF (MODI N ITER,NPERR).EQ.O)  WRITE<6,202)  NITER, SUM 
202  FORMAT ( 20H0  NO.  OF  I T6RATI ONS=I 4 , 5X,6HERR0R=E15 .7 > 

I F ( MOO (NITER,NPOTPT).EQ.O)  CALL  OUTD I S ( U , H, I MAX , JMAX ,KMAX , BO IM, 
1  CDIM, 1 ) 

C  CHECK  IF  ERROR  MEETS  TOLERANCE 

IF  (SUM-TOL)  12,12,13 
13  IFINITER.LT. NIT)  GO  TO  10 

12  WR I TE ( 6, 202 )  NITER, SUM 
NUPRNT=NUPRNTG1 

GO  TO  (401,890,900) ,NUPRNT 

890  WRITE  (10)  ( (  ( U(  I ,  J  »K)-»  K=1 » 30 )  ,  J— -l-»-10-)~,-I*  1 ,  IMAX  ) 

WRITE  (6,891) 

891  FORMAT  ( 1HO,5X,30HDI SPLACEMENTS  SAVED  ON  TAPE  10  > 


TPA00920 

TPA00930 

TPA00960 

TPA00970 

TPA00980 

TPA00990 

TPA01000 


TPA01020 


TPA  103 

TPA01040 

TPA01050 

TPA01060 

TPA01070 

TPA01080 

TPA01090 

TPA01100 

TPA01 1 10 

TPA01120 

TPA01130 

TPA01140 

TPA01150 

TPA01160 

TPA01170 

TPA01180 

TPA01190 

TPA01200 

TPA01210 

TPA01220 

TPA01230 

TPA01240 

TPA01250 

TPA01260 

TPA01280 

TPA01290 

TPA01310 

TPA01320 
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co 


GO  TO  401 

900  READ  (5,902)  I , IMAX  ,  JMAX ,KMAX ,N I T ,NPERR, NPOTPT ,T0L ,XFAC 
902  FORMAT  (7I5.2E10.3) 

IF  ( I.NE.(-619) )  GO  TO  401 
CALL  INTERP 
GO  TO  300 

CALCULATE  AND  OUTPUT  STRESSES 
401  CALL  PART2 
GO  TO  1 
END, 

SUBROUTINE  0 1 SPL (H , T, ALF A ,0 1 MJ ,0 1 MK , I M, JM,KM,U, D ) 


INITIAL  ESTIMATE  OF  DISPLACEMENTS 


TPA01330 

TPA01340 

TPA01350 

TPA01360 

TPA01370 

TPA01410 

TPA01420 

TPA01430 

TPA01440 

TPA01450 

TPA01460 

TPA01510 

TPA01480 

TPA01490 

TPA01500 


DIMENSION  H(20) ,T(20) , ALFA (20) ,U (21 , 10, 30 > , XI ( 21 ) , Z J (2 1 ) , YK ( 21 ) , V ( 
110),Z(10) 


NL* IM— 1 
DO  1.  1*1, NL 

2 J ( I )  =  T  (  I ) * ALFA ( I  ) *DI M J 

1  YK( I )*T( I)*ALFA( I )*DIMK 
XI(1)=0 

DO  2  1*2, IM 

XIII ) -x I  ( I-DCHl  I-1)*ALFA<  I-1)*T(I-1) 
ZJ(I)*(ZJ(I-l)£ZJ(I)>/2. 

2  YK(I)»(YK(I-l)£YK(I))/2i 
YK(IM)=T(NL)»ALFA(NL)*DIMK 
ZJ ( IM )=T (NL )* ALFA (NL ) *DI MJ 
Y(i>=0. 

2(1) -0- 
DO  3  1=2, JM 

3  Z(I)=FLOAT(I-l)/FL0AT(JM-l) 

DO  4  1=2, KM 

4  Y( I )=FLQAT ( I— 1 ) /FLOAT (KM— 1 ) 

DO  5  1  =  1 ,  IM 

DO  5  J*  1  * JM 
00  5  K=1 » KM 
U( I , J ,K ) =XI ( I ) 

U( I, J,KC10)=YK( I )*Y(K) 

5  U( I , J ,KG20) =2 J ( I ) * Z  ( J ) 


TPA01540 

TPA01550 

TPA01560 

TPA01570 

TPA01580 

TPA01590 

TPA01600 

TPA01610 

TPA01620 

TPA01630 

TPA01640 

TPA01650 

TPA01660 

TPA01670 

TPA01680 

TPA01690 

TPA01700 

TPA01710 

TPA01720 

TPA01730 

TP A  175 


S=0. 

DO  6  1=1, IM 
DK=U( I, JM, KMC 10) /2.0 
IF  ( 1-1 )  61,61,62 
61  D$*H( I )*DK 


TPA01770 

TPA01780 

TPA01800 

TPA01810 
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non 


GO  TO  6 

TPA01820 

62 

DS*( HI  1 ) CHI 1-1) )*DK 

TPA01830 

IFII.EQ.IM)DS=HI I - 1 ) *DK 

TPA01840 

6 

S-SSDS 

TPA01850 

DEL* ( S/D-UJ IM.JM.KMS10) )/D*DIMK 

TPA  186 

00  8  J-l.JM 

TPA01870 

Z-FLOAT ( J- I ) /FLOAT I JM- 1 ) *DI M J 

TPA01880 

DO  8  K*  1  »KM 

TPA01890 

Y®  FLOAT ( K-l )/ FLOAT ( KM-1 ) *DI  MK 

TPA01900 

R*(Z«Z£Y*Y)/D1MK*»2*DEL 

TPA01910 

DO  8  1*1, IM 

TPA01920 

8 

U<  I,J,K)*UI I, J»K)£R 

IMlD*IM/2 

UM  ID*U ( IMIOfl.l) 

TPA01940 

DO  9  1*1 » IM 

TPA01960 

DO  9  J-l.JM 

TPA01970 

DO  9  K*1,KM 

TPA01980 

9 

U(  I»J,K)=U( I  *  J t  K )-UM ID 

RETURN 

TPA02000 

END 

TPA02010 

SUBROUTINE  SOLVE (A.L.M.N) 

TPA02060 

#*##1 

GAUSS  SEIDEL  ITERATION 

TPA02040 

DIMENSION  A(18.9f9) , FR ( 3 ) ,UR ( 3 ) , D ( 3 ) • DU ( 3 ) 

COMMON/SOL/F( 21, 10,30} ,U (2 1 , 10, 30) , I  MAX , JMAX.KMAX 

, SUM, XF AC, TEMPI  20 

l ) »  ALFA ( 20 ) » H  ( 20) , E ( 20 ) , PR ( 20 )  ,BDIM,CDIM 

COMMON/ 1 DEN /Ml DENI 10) ,NIDEN<10> .NITER 

DOUBLE  PRECISION  F 

DOUBLE  PRECISION  FR 

TPA02110 

IK*3*MIDENIM)-3£NI0ENIN> 

TPA02140 

DO  1  1*1,3 

TPA02150 

I 3=N£  l I- 1 ) • 10 

I4=I£t I  —  1 ) »  3 

URI1 )*UIL.M,I3) 

FR ( I ) *F ( L ,M , 1 3) 

TPA  216 

I 

0( I )=A( IK, IK, 14) 

MML*M- 1 

TPA02L90 

MMU*M£1 

TPA02200 

NNL*N- 1 

TPA02210 

NNU=N£1 

TPA02220 

IF  IMIDEN(M)-2)  5,7,6 

TPA02230 

5 

MML*1 

TPA02240 

GO  TO  7  - 

TPA02250 

6 

MMU® UMAX 

TPA02260 
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non 


7  IF  (NIDEN(N)-2)  8,10,9 

8  NNL=1 

GO  TO  10 

9  NNU=KMAX 
10  NNLTMP=N 

MMLTMP=M 

LL=L 

19  DO  21  MM=MMLTMP,MMU 
DO  20  NN=NNITMP,NNU 

IJ  =  3» (MI DEN (M) CMM-M) GNIDEN(N) GNN-N-3 

15  DO  16  I X=1 , 3 

DO  16  JX*1,3 

1 3*IX£( JX-1 ) *  3 
I4*NN£(JX-1)*10 

16  FR(IX)=FR(I X) —A (IK,IJ,I3) *U ( LL  »MM *  14 ) 

20  CONTINUE 

21  NNLTMP*NNL 

IF  < IK.GT.9)  GO  TO  24 

IF  (L.EQ.IMAX)  GO  TO  24 

LL*LLG1 

IK*IK.£9 

MMLTMP=MML 

m  GO  TO  19 

24  CONTINUE 
DO  70  I 

50  FR ( I )*FR (I ) CD ( I ) *UR ( I ) 

0U( I ) =FR  Cl) /D II) -UR ( I ) 

70  DU( I >=DU( I ) *XFAC 
IF(N.EQ.l)  DU ( 2 ) =0 
IF(H.EQ.l)  OU ( 3 ) =0 
DO  80  1=1,3 

SUM=SUMCABS (DU( I ) *D ( I ) /XFAC ) 

I3=NC( 1-11*10 

80  U(L,M,I3)SU(L,M,I3) CDU ( I ) 


MODIFY  FORCES  TO  ACCOUNT  FOR  SYMMETRY 


DO  27  1*1,3 
I3*Nfi( 1—11*10 

27  DU( I ) =U( L ,M, I  3 ) 

28  LL=l 

IK*3«MIDEN(M)-3CNIDEN(N) 

MHLTMP*M 

NNLTHP=N£1 

IF  (NNLTMP.LE.NNU)  GO  TO  29 


TPA02270 
TPA02280 
TPA02290 
TPA02300 
TPA02310 
TP ‘02320 
TPA02330 
TPA02340 
TPA02350 
TPA02360 
TPA02370 
TPA02380 

TPA  240 

TPA02400 
TPA02410 
TPA02420 
TPA02430 
TPA02440 
TPA02450 
TPA02460 
TPA02470 
TPA02480 
TPA02490 
TPA02500 
TPA025 10 
TPA02520 
TPA02530 
TPA02540 
TPA02550 
TPA02560 
TPA  257 

TPA02580 
TPA02590 
TPA02600 
TPA02610 
TPA  262 

TPA02630 

TPA02640 

TPA02650 

TPA02660 

TPA02670 
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c 

c 


NNLTMP=NNL 
MML  TMP=M£1 

IF  (MMLTMP.GT.MMU)  GO  TO  31 

29  DO  30  MM=MMLTMP,MMU 
DO  26  NN=NNLTMP , NNU 

IJ  =  3» ( MID5N (M) CMM-M ) CNIDEN! N) 6NN-N— 3 

25  00  26  IX=1,3 

DO  26  JX-1 *3 
I4«JX£( I X-l ) »  3 
I3=NNG( IX-1)*10 

26  F(LL,MM, 13) *F ( LL ,MM, 13 ) -A ( IK, I J,I4)*DU(JX) 

30  NNLTMP=NNL 

IF  ( IK.GT.9)  GO  TO  35 

31  IF  (L.EQ.IMAX)  GO  TO  35 
LL*LL6l 

IK* IKG9 
MMLTMP-MML 
GO  TO  29 
35  CONTINUE 
900  RETURN 
END 

SUBROUTINE  OUTDIS  ( U,H , I  MAX » JMAX ,KMAX »BDI M , CDI M, NT1 ) 

PRINTS  OUT  DISPLACEMENTS  AND  COORDINATES  OF  NODES 
DIMENSION  XC<  3) ,U(21 .10*30) »H(20> 

K0UNT=3 

GO  TO  ( 10f 20) »NT1 
10  WR I TE ( 6 1 900 ) 

900  FORMAT! 1H1.35X.29HD  ISPLACEMENTS  U// 10X  *  1HX  » 1 IX  * 1HY* 
l  11X, 1HZ, 16X,3HU-X,14X,3HU-Y» 14X.3HU-Z) 

GO  TO  30 
20  WRITE  (6,899) 

899  FORMAT! 1H1.35X, 16HF  0  R  C  E  S  F  //10X, 1HX, 11X, 1HY, 

1  11X, 1HZ,16X,3HF-X, 14X , 3HF-Y, 14X , 3HF-Z ) 

30  XC( 1 )=0. 

00  101  L=1 » IMAX 
WR I TE ( 6, 901 ) 

K0UNT=K0UNTC1 
DO  100  K=1 ,KMAX 
XC(2)=FL0AT(K-1)*BDIM 
DO  100  J=1 , JMAX 
XC(3)=FL0AT(J-1)*CDIM 

WRITE (6, 902)  XC,U(L,J,K) ,U(L,J,KC10) ,U(L, J ,KC20 ) 

K0UNT*K0UNTE1 
IFIKOUNT.lt. 51)  GO  TO  100 


TPA02680 

TPA02690 

TPA02700 

TPA02710 

TPA02720 

TPA02730 

TPA02740 

TPA02750 

TPA  277 

TPA02770 

TPA02780 

TPA02790 

TPA02800 

TPA02810 

TPA02620 

TPA02830 

TPA02840 

TPA02850 

TPA02660 

TPA03050 

TPA03060 

TPA03070 

TPA03090 

TPA03100 

TPA031 10 

TPA03120 

TPA03130 

TPA03140 

TPA03150 

TPA03160 

TPA03170 

TPA03180 

TPA03190 

TPA03200 

TPA03210 

TPA03220 

TPA03230 

TPA03240 

TPA03250 

TPA  402 

TPA03270 

TPA03280 
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K0UNT=0 

GO  TO  (  110,  120)  » NT1 
110  WR I TE ( 6 » 900 ) 

GO  TO  100 
120  WRITE  ( 6*899) 

100  CONTINUE 

101  XC(1)=XC<1)£H(L) 

RETURN 

901  FORMAT ( 1H0) 

902  FORMAT! IX, 3F1 2.6 , 5X , 3 ( 1PE1 7. 5 ) ) 
END 

SUBROUTINE  ASSEM(A,AT,BIK) 


ASSEMBLE  STIFFNESS  MATRIX  BY  LAYERS 


DIMENSION  A(18,9,9) ,AT(18,9,9) 
DIMENSION  I LOC (16) 

EQUIVALENCE  ( I LOC ( 1 ) , LOC 1 1 , 1 ) ) 

DATA  (ILOC(J) »J=l,16)/4*5*7,8» 

DO  1  1*1, A 
DO  1  J=  1  *  4 
JB-LOCC I,J) 

JT*JB£9  . 

JJ*J£4 
DO  1  K*1 ,4 
J  5* J J  £ ( K- 1 ) *8 
J4=J£ ( K— 1 ) *  8 
KB=LOC ( I , K ) 

KK*K£4 

J3=JJ£( KK-1  )*8 
DO  1  L®1 ,3 
DO  1  M*1 , 3 
I3*L£(M-1)*3 
A ( JB ,KB , 1 3 ) *A ( JB ,KB , 13 ) £BLK ( L ,M, J3 ) 

A I JT , KB , 1 3 ) =A ( J  T , KB , 1 3 ) CBLK ( L , M , J5 ) 

1  AT(JB,KB,I3)=AT(JB,KB,I3)£BLK(L,M,J4) 
RETURN 
END 

SUBROUTINE  FSEQ 


FORM  SYSTEM  OF  EQUATIONS 


, LOC (4, 4) ,BLK(3,3,64) 
5,6,8,9,2,3,5,6,1,2,4,57 


TPA03290 
TPA03300 
TPA03310 
TPA03  320 
TPA03330 
TPA03340 
TPA03350 
TPA03360 
TPA03370 
TPA03380 
TPA03390 
TPA03700 
TPA03670 
TPA036B0 
TPA03690 


TPA03730 

TPA03740 

TPA03750 

TPA03760 

TPA03770 

TPA03780 


TPA03790 

TPA03800 

TPA03810 

TPA03820 


TPA03860 

TPA03870 

TPA03920 

TPA03890 

TPA03900 

TPA03910 


COMMON  AB ( 18, 9, 189 ) 

COMMON/SOL 7F( 21, 10,30) ,U(21 , 10,30) , IMAX, JMAX,KMAX,SUM,XFAC,TEMP( 20 
1), ALFA (20) ,H(20) , E ( 20) , PR ( 20 ) ,BDIM,CDIM 
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o  o  o 


o 

to 


DOUBLE  PRECISION  F 

DIMENSION  8LK (3.3.64)  ,BLKT 1 3 ,3 ,64 J , AT (18 *9,9) 

EQUIVALENCE  ( U  ( l ,  I ,  I )  ,  BLM  1 , 1 1 1 ))  .  ( U(  1 1 1 »  10  )  fBLKT  ( 1 , 1,  I  >) 

DO  I  1=1,21 

DO  I  J=l,10 

DO  I  K  =  1 » 10 

DO  1  L-l , 3 

I3*K£(L-1)*10 

F ( I i J  « I  3 ) *0. 

1  UII,J,I3)=0. 

00  2  1-1,18 
DO  2  J  =  l,9 
DO  2  K=1 , 3 
DO  2  L*i  ,3 
I3»K£(L-i)*3 
ATI  I , J , 1 3 ) =0. 

DO  2  M=1 » IMAX 

I 4=K£I L— 1 )*3£ I M- l ) *9 

2  AB( I, J, I4)=0. 

LMAX= IMAX— 1 
DO  3  L-l ,LMAX 
Ll=L£l 

ALFAT=TEMPIL)*ALFAIL) 

CALL  BSTIFIH(L) ,BDIM ,CDI M,E ( L ) , PR { L ) ,BLK , BLKT ) 

CALL  ASSEMI AB l 1 , 1 »9*L-B ) , AT ,BLK ) 

CALL  LOAD(H(L),BDIMtCDIM,E(L)fPR(L) ,ALFAT, JMAX,KMAX,L,F ) 
DO  4  1=1,18 
00  4  J*l»9 
DO  4  K=l, 3 
DO  4  LK=1,3 

I3*K£( LK-1 )*3£(Ll-l)*9 
I4-KS ( LK- 1 )*3 

AB( I , J, 13) *AB 1 1 ,J, 1 3) CAT (I , J, 14) 

4  ATI I,J,I4)»0. 

3  CONTINUE 
RETURN 
END 

SUBROUTINE  BFORCE (A,B,C,E,PR,ALFAT»FB) 


COMPUTE  THERMAL  NODE  FORCES  ON  BLOCK 


DIMENSION  FB(3,B) ?BDIM(3),IPHI (8,4) ,IPHITE(32) 
EQUIVALENCE  (IPHlTE(i) , I  PH I ( 1 • 1 ) ) 

DATA  I IPHITEI J), J=l,32)/  Cl , £1 ,B1 , Cl ,-L ,-l ,-l ,-l , 

l  -l,£lt£lf-l,-l,£l,£l,-l. 


TPA04010 


TPA04040 

TPA04050 

TPA04060 

TPA04070 


TPA04090 


TPA041 10 
TPA04120 
TPA04130 
TPA04140 
TPA04150 

TPA04170 

TPA04190 

TPA04200 

TPA04210 

TPA04220 


TPA04250 

TPA04260 

TPA04270 

TPA04320 

TPA04290 

TPA04300 

TPA04310 

TPA04330 
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non 


2 

3 


£l,£i»-l»-l»El,£l»-l,— 1, 
-1, £1,-1, El, £1,-1, £1,-1  / 


o 

w 


Cl*£/< l.-2.«PR> 

BDIH ( 1 ) 3  A 
BOIMI 2)=B 
BD1M(3)-C 
V=A*B*C 
00  1  1-1,3 
CI-V/4./BDIMCI  ) 

00  1  J-1,8 
IS«1 

00  2  M=l, 4 
IF(M.EQ.I)  GO  TO  2 
IS*IS*IPHI ( J, M) 

2  CONTINUE 
S-IS 

1  FBI I,J)*C1*CI*ALFAT*$ 

RETURN 

END 

SUBROUTINE  LOAD ( A, B  ,C ,E  ,PR , ALFAT •  JHAX « KMAX , L, F ) 


ADD  THERMAL  FORCES  TO  NODES  IN  LAYER 


DIMENSION  F(21fl0,30),FB(3y8), LOC ( 8 ) 

DOUBLE  PRECISION  F 

DATA! LOC (I ), 1-1,81/8,7,5,6,4,3,1, 2/ 

JU-JMAX— 1 

KU*KMAX— 1 

LL*L 

LU-LE1 

CALL  BFORCEIA,B,C,E, PR, ALFAT, FBI 
DO  1  J-i,JU 
DO  1  K=1 ,KU 
KK-0 

00  1  IL=LL,LU 
DO  1  M-1,2 
00  l  N-l»2 
KK-KKE1 
JF-LOC(KK) 

JM-JEM-1 
KN-KEN-1 
DO  1  IN-1,3 
I 3=KNE ( IN- 1 1* 10 

1  F(IL,JM,I3)*F( IL* JM, 1 3) EFB ( IN*  JF ) 
RETURN 


TPA04360 

TPA04370 

TPA04380 

TPA04390 

TPA04400 

TPA04410 

TPA04420 

TPA04430 

TPA04440 

TPA04450 

TPA04460 

TPA04470 

TPA04480 

TPA04490 

TPA04500 

TPA04510 

TPA04520 

TPA04570 

TPA04540 

TPA04550 

TPA04560 


TPA04590 
TPA04600 
TPA04610 
TPA04620 
TPA04630 
TPA04640 
TPA04650 
TPA04660 
TPA04670 
TPA04680 
TPA04690 
TPA04700 
TPA04710 
TPA04720 
TPA04730 
TPA04740 
TPA04750 
TPA  476 

TPA04770 
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ENO 

SUBROUTINE  INTERP 

C  INTERPOLATES  DISPLACEMENTS  U  TO  A  GRID  OP  DOUBLE  SIZE 

C  RESULT  IS  INSERTED  INTO  F 

COMMON/ SOL/FI  21, 10,30) ,U(21 ,10,30) , I  MAX, JMAX, KMAX , SUM.XFAC, TEMP ( 
1), ALFA (20), HI  20) ,E ( 20) ,PR(20) ,BDIM,CDIM 
DOUBLE  PRECISION  F  -  ' 

DO  10  L= 1,1 MAX 
00  10  J-l , JMAX 
Jl*(  JED/2 
J2  =  l 

IF  (  (MOD  (J,2)).EQ.O 


TPA04780 

TPA05350 

TPA05360 

TPA05370 


20 


)  J2=2 


o 


10 


J3=J  1£ J2-1 
DO  10  K=1,KMAX 
Kl*(K6l)/2 
K2=l 

IF  (  ( MOD ( K , 2 ) ) . EQ . 0  )  K2=2 
K3=K1CK2-1 
DO  6  1*1,3 

S  =  0. 

DO  5  JZ=J1? J3 

DO  5  KZ=Kl,K3 

I3=KZ£( I- 1 ) *  1 0 

1  4=K  E  ( I— 1 ) *  10 

S*SCU(L, JZ, 13) 

F(L,J,IA)*S/FLOAT(K2*J2) 

IF(K.EQ.l)  F(L,J,KE10)*0. 

IF  (J.EQ.l)  F (L, J,K£20) =0. 

CONTINUE 

RETURN 

END 

SUBROUTINE  PART2 


THERMAL  STRESS  ANALYSIS  OF 
PART  II  STRESS  OUTPUT 


RECTANGULAR  PRISM 


TPA05400 

TPA05410 

TPA05420 

TPA05430 

TPA05440 

TPA05450 

TPA05460 

TPA05470 

TPA05480 

TPA05490 

TPA05500 

TPA05510 

TPA05520 

TPA05530 

TPA05540 


COMMON/STR/  XC1 ,KOUNT, J ,K 

COMMON/SOL/F( 21, 10,30) ,U( 2 1 , 10 , 30 ) , I  MAX , JMAX ,KMAX ,SUM,XFAC, TEMPI 
1),ALFA(20),H(20) ,E(20),PR(20) ,BDIM,CDIM 
DOUBLE  PRECISION  F 

COMMON/ I  DEN/M  I  DENI  10) ,NIDEN( 10) »NI TER 
DIMENSION  UT( 8 , 3 ) 

PRINTOUT  DISPLACEMENTS 


TPA05590 

TPA05600 

TPA05610 

TPA05640 

,TPA05650 

TPA05660 

TPA05670 

»TPA05680 

TPA05690 


20 


TPA05730 

TPA05750 

TPA05760 

TPA05770 
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CALL  OUTOIS  (U.H.IMAX, JMAX , KMAX , BDI M ,CDI M, 1 ) 

PRINTOUT  STRESSES  AT  CENTR0I0  OF  EACH  ELEMENT 


WR I TE  (  6 . 1000) 

1000  FORMAT I 1H1.30X.65HS  T 
1  ELEMENT  S// ) 
KOUNT=3 
I M= I MAX- 1 
JM=JMAX- 1 
KM=KMAX- I 
XCl=H(l)/2. 

DO  155  L=1 » IM 
ALFAT=ALFAI L ) *TEMP  (  L ) 
DO  150  K=1 « KM 
DO  150  J*1 . JM 
CALL  S£TUT(UT,L, J.K) 


R  E  S  S  E  S 


A  T 


CENTROIDS 


0  F 


CALL  j.lul  HI!  .  L  .  J  «  A  I 

150  CALL  STRES(H(L).8DIM.CDIM.0.«0.«0a.UT.E(L).PR(L).ALFAT) 
155  XC1=XC16IH(L61)6HIL) 1/2.0 


900  RETURN 
END 

SUBROUTINE 


STRES  (A.B.C.PSI .ETA.ZETA.U.E.PR.ALFAT ) 


SUPPLY  U(8, 31-DISPLACEMENTS  AT  8  NODES  OF  ELEMENT 

A,B,C  -  DIMENSIONS  OF  ELEMENT  E.PR-  ELAST  CONSTANTS 
PSI »ETA»  ZETA  -  LOCAL  COORD.  WHERE  STRESS  IS  TO  BE  CALCD 

S ( 7 )  AND  EPSII7)  WILL  GIVE  STRESSES  AND  STRAINS  RESP. 

(XX.YY.ZZ.XY.XZ.YZ,  AND  EFFECTIVE ) 

5  AND  PRINCIPAL  STRESSES  AND  DIRECTIONS  WILL  BE  PRINTED  OUT 
KOUNT  IS  LINE  COUNT  ON  PAGE  BEING  PRINTED 
XC1  IS  X-COORDINATE  OF  CENTROID  OF  LAYER 

COMMON/ STR /  XC 1 , KOUNT , JZ , KZ 
DIMENSION  XC ( 3 ) , W1 ( 6 ) , W2 ( 3 , 3 ) ,W3(3> 

DIMENSION  FEE  (8,4)  ,D  ( 3 1  ,  XLOC  I  3  )  ,U  18 ,3 )  .EPS 1 1 7  »  »-S(7)  ,FEEID<  32) 
EQUIVALENCE  C FEE  ID  Cl) .FEE C 1  ,1 1) 

DATA (FEE  ID ( J) , J* 1, 32) /El. ,C1. ,GI. ,GI. .-1..-1..-1 .,-1.. 

1  &l.»£l.»-l.»-l.»£l.»tl.»-l.»-I.» 

2  — l.»EI.»El«»  — 1»»  — —  1.. 

3 


TPA05780 
TPA05790 
TPA05800 
TPA05810 
TPA05820 
TPA05830 
TPA058A0 
TPA05850 
TPA05860 
TPA05870 
TPA05880 
TPA05890 
TPA05900 
TPA05910 
TPA05920 
TPA05930 
TPA05940 
TPA05950 
TPA05960 
TPA062 1 0 
TPA06220 
TPA06240 
TPA06250 
TPA06260 
TPA06270 
•TPA06280 
TPA06290 
TPA06300 
TPA06310 
TPA06320 
TPA06330 
TPA06340 
TPA06350 
TPA06360 


D( 1 ) =1 . /A 
D(2I=I./B 
D(3)=l./C 
XLOC ( 1 ) sPS I 


TPA06420 

TPA06430 

TPA06440 

TPA06450 
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XLOC ( 2) -ETA 
XLOC ( 3 ) =ZETA 
C 

1 STRAN=1 
00  21  11*1,3 
00  21  12*1,11 
SUM1=0. 

SUM2=0. 

00  20  J= 1 » 8 
PR001=FEE( J,4» 

PR0D2*FEEf J.4) 

00  12  K*l,3 
IF  (K.EQ.I2)  GO  TO  11 
PROD1=PROD1*(XLOCIK)GFEE(J,K) ) 

11  IF  (K.EQ.I1)  GO  TO  12 
PROD2*PROD2*(XLOC(K) GFEE ( J,K) ) 

12  CONTINUE 

SUM1-SUM 1 CPROO 1 *U ( J  » 1 1 ) *0 ( 1 2 ) 
SUM2=SUM2£PROD2*U( J • 1 2 1 *0 (III 

20  CONTINUE 

S  < ISTRAN>=(SUM16SUM2)/8.0 

21  ISTRAN*ISTRANC1 

Epsim*sm 

EPSI ( 2)=S( 3) 

EPSI ( 3) =S( 6) 

EPSI ( 4) *Sl 2 ) 

EPSI ( 5 )=S( 4 ) 

EPSI ( 6  >  =S ( 5  > 

00  25  1=1.3 

25  EPSIII 1*EPSI ( I J-ALFAT 

DIL  =  EPS I ( 1 ) CEPS  I (2) CEP SI (3) 
Cl*E*PR/(  l.£PR)/U.-2.*PR> 
C2*E/(1.CPR> 

SUM2=0. 

00  30  J=l,3 
S(J)=C1*DIL£C2*EPSI (J) 

30  SUM2=$UM2CSIJ> 

00  31  J =4  *  6 

31  S(J»=C2*EPSI( J) 

C  CALCULATE  EFFECTIVE  STRESS  AND 

SUM1=0. 

00  40  J*l»3 

40  SUM1=SUM16IS( JJ-SUM2/3. )*»2 
00  41  J  *4 . 6 

41  SUM1*SUM162.*S( J>**2 


TPA06460 
TPA06470 
TPA06480 
TPA06490 
TPA06500 
TPA06510 
TPA06520 
TPA06530 
TPA06540 
TPA06550 
TPA06560 
TPA06570 
TPA06580 
TPA06590 
TPA06600 
TPA06610 
TPA06620 
TPA06630 
TPA06640 
TPA06650 
TPA06660 
TPA06670 
TPA06680 
TPA06690 
TPA06700 
TPA06710 
TPA06720 
TPA06730 
TPA06740 
TPA06750 
TPA06760 
TPA06770 
TPA06780 
TPA06790 
TPA06800 
TPA06810 
f PA06820 
TPA06830 
TPA06840 

STRAIN  TPA06850 

TPA06860 

TPA06870 

TPA06880 

TPA06890 

TPA06900 
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c 

c 


S I 7 ) =SQRT( SUM  1 ) 

TPA06910 

TPA06920 

OBTAIN  PRINCIPAL  STRESSES  AND  DIRECTIONS 

TPA06930 

Ml ( 1)=S( 1) 

Mil  2 ) =S I 4 ) 

Wl( 3 )=S( 2) 

W H  4  J  =S ( 5) 

W1 ( 5 1 =S( 6 ) 

W1 ( 6 ) =S( 3 ) 

CALL  JACOBI (W1,W2,W3, 3) 

IF  ( KQUNT.GT . 5 )  GO  TO  60 

WR I TE ( 6 , 101 ) 

TPA07090 

101  FORMAT! 2 5X, 74HP0SI T ION  /  STRESS  COMPONENTS  /  PRINCIPAL  STRESSES  ANTPA07100 
ID  DIRECTIONS  «3  L I NES )  /  14X ,  1HX  ,  1 3  X,  1HY ,  1 3X ,  1HZ/  12X  ♦  2HXX  ,'15X  ,  2HYY  »  TPA07110 
215X,2HZZ,15X,2HXY»15X,2HXZ,15X,2HYZ,9X»10HEFF. STRESS/25  X(,  TPA07120 

2  12HPRINC. STRESS, 11X.  6HX-COMTPA07 130 


o 

-a 


3P » 1 IX, 6HY-C0MP , 1 1 X 1 6HZ-C0MP  * 3X  *  1 3H0F  PR  INC. DIR. ) 
KOUNT=KOUNT£3 
60  XC< 1  )=XCl£PSI*A/2. 

XC(2)=FL0AT(KZ-l)*B£(ETA£1.0)*B/2. 

XC(3)=FLOAT(JZ-i)*C£(ZETA£1.0)*C/2.0 

WR I TE ( 6 , 102  )  XC,S,(I,W3(IJ  , I W2 ( J , II , J  = 1 ,3 ) , I* 1 , 3 ) 

102  FORMAT ( 1HO,  3X  ,  3F  14 . 6/ 1 X ,  7  ( l PE  1 7 . 5  >  / 1 15X , 1 6 » 1PE 17  .5 f 3F17 . 5  )  ) 
K0UNT=K0UNT£6 

IFIKOUNT.LT. 53)  GO  TO  70 

K0UNT=0 

WRITEI6, 103) 

103  FORMAT  ( 1H1) 

70  CONTINUE 

RETURN 

END 

SUBROUTINE 


SETUT ( UT ( L • J • K ) 


TPA07140 

TPA07150 

TPA07160 

TPA07170 

TPA07180 

TPA07190 

TPA07200 

TPA07210 

TPA07220 

TPA07230 

TPA07240 

TPA07250 

TPA07260 

TPA07270 

TPA07280 

TPA07300 

TPA07310 

TPA07320 

TPA07330 


OBTAINS  DISPLACEMENTS  UT  FOR  A  PARTICULAR  ELEMENT 

COMMON/ SOL/FI  21, 10,30) »UI 21 » 10,30) « I MAX * JMAX , KMAX « SUM, XF AC, TEMP  I  20 
1) , ALFA (20) ,HI 20) ,E(20) ,PR(20) ,BDIM,CDIM 
DOUBLE  PRECISION  F 

DIMENSION  UT( 8* 3 )  TPA07360 

L2*L£1  TPA07370 

L3*l  TPA07380 

DO  40  L 1=1 • 2  TPA07390 

DO  30  1=1,3  TPA07400 

I3*KCl I- 1 >  *  10 
UT(L3, I ) =U ( L2 , J , I 3G1 ) 
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noooo  non  ooonoonno 


UT(L3C1,I)=U(L2,J61,I3£1) 

UT(L3C2,I)=U(L2,J£l,I3) 

30  UT(L3£3, I)=U(L2, J, 131 
L2=L2-1 
40  L3=L3C4 
RETURN 
END 

SUBROUTINE  JAC08I (A,B,C»N) 
DIMENSION  A ( 3 ) t B ( 3 «  3 ) >C(3) 


TPA07450 

TPA07460 

TPA07470 

TPA07480 


N  =  ORDER  OF  MATRIX 
A  *  INPUT  MATRIX 

VECTORS  FORMED  AND  RETURNED  IN  ARRAY  B 
VALUES  RETURNEO  IN  C 


SET  UP  IDENTITY  MATRIX  FOR  VECTORS 

00  15  J=1,N 
00  15  1=1, N 
IF(I-J)  10,5,10 
5  B ( I , J ) =1 .0 
GO  TO  15 
10  B( I , J ) =0.0  - 
15  CONTINUE 

CALCULATE  NORM  OF  MATRIX 

XNORM=0.0 

IJ=0 

DO  30  J=2,N 
I  J  =  I  J£1 
JM1=J-1 
00  30  1=1, JMl 
IJ=IJ£1 

30  XNORM=XNORMEA ( I JI*A(I J) 

XN0RM=SQRT (2. 0*XN0RM) 

IF(XNORM)  120,120,32 

SET  CONVERGENCE  CRITERION 

CONVERGENCE  WHEN  LARGEST  OFF-DIAGONAL  ELEMENT  IS  LESS  THAN 
2**-27  TIMES  NORM  OF  INPUT  MATRIX 

32  ZNORM=XNORM*< . 74505806E-08) 

SIGMA=  FL0AT(N)*1. 414214 
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nno  ooo  non  n  n  o  onnnnon 


LOWER  THRESHHOLO  BY  2»*l/2  TIMES  N  WHEN  NO  OFF-DIAGONAL 
ELEMENTS  EXCEED  CURRENT  VALUE 

35  XNORM=XNORM/S IGMA 

NEW  SWEEP  THROUGH  MATRIX 

AO  IN0=0 
IJ*0 

A5  00  100  J  =  2*N 
IJ»IJ£1 
JM 1= J- 1 

DO  100  1*1,  JM 1 
I J=IJE1 

IF(ABS(A( I J) J-XNORM)  100,100,50 

ELEMENT  EXCEEDS  THRESHHOLO  -  NEW  SWEEP  NECESSARY 
50  IND=1 

PIVOTAL  SET  FOR  THIS  ROTATION 

II  =  < I*(I£1) )/2 
JJ-( J*(J£1) )/2 
TEMP=A(IJ) 

0 1 1  -  A  ( 1 1  ) 

DJJ=A( J J ) 

XLAM=-TEMP 

XMU=0.5*(DII-DJJ) 

XMEGA*. XLAM/SQRT(XLAM*XLAM£XMU*XMU) 

IF(XMU)  55,60,60 
55  XMEGA=-XMEGA 

CALCULATE  SINE  AND  COSINE  OF  ANGLE  OF  ROTATION 

60  SINE*  XMEGA/I SQRT(2 . 0* ( 1 . OGSORT ( 1 . 0-XMEGA*XMEGA) ) )) 
COSIN=  SORT (1.0-(SINE**2)) 

TRANSFORM  MATRIX  AND  VECTORS 

kj*(  j*i  j-m/2 
K 1  =  1 1* ( I- 1 ) ) /2 
00  95  K-lt N 
IF(K-I)  85,62,65 
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non 


62  K 1  =  I  I 
KJ-K J£1 
GO  TO  92 

65  IF(K-J)  70,67,80 
67  KJ«JJ 

K I =K I £K—  1 
GO  TU  92 
70  KJ=KJ£1 
75  K I =K IfiK-  1 
GO  TO  90 
80  KJ=KJEK-L 
GO  TO  75 
85  KJ=KJ£1 
KI=KI£1 
90  AK  J=A ( K J ) 

A ( K J ) *  A(KI )*SINE£AKJ«C0SIN 
A ( K I  )  —  A ( K I )*C0$I N-AK J*S I NE 
92  BKI*  B ( K  ,  I ) 

B(K,I)  =  8K I *COS I N-B ( K , J ) *S I NE 
B(K,JJ*  BKI*SINE£B(K»J)*C0SIN 
95  CONTINUE 

A( I J )=  (DII-DJJ)*(CQSIN»SINE)£T£MP*( (COS IN*»2)-(SINE**2)> 

A( II )=  DII*(COSIN**2)£0JJ*(SINE»»2)-(2.0*TEMP*<C0SIN*SINE) > 
A ( J J ) =  OII*(SINE«*2)£OJJ*(COSIN»*2)£(2.0*TEMP»(COSIN*SINE) ) 
100  CONTINUE 

IF(IND)  AO, IIO, AO 
110  IF ( XNORM-ZNORM)  120,35,35 


SORT  VALUES  AND  VECTORS 


120  00  1 AO  1=1, N 
II*( I«( I £1 )  )/2 
TEMP=A( II) 

00  130  JSI «N 
JJ»( J*( J£l)  )/2 
IF ( TEMP-A ( J J )  )  130,125,125 
125  IT*J 
ITJ*JJ 
TEMP-A (JJ) 

130  CONTINUE 

C( I )*A( ITJ  ) 

A(ITJ)»A(I  I ) 

00  1A0  K=1  »N 
A( II)*B(K,I ) 

B(K, I )=B(K, IT) 
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non 


140  8(K« IT) = A (  II) 


NORMALIZE  VECTORS 

DO  L75  J=ifN 
XNORM^O.O 
00  155  I =1 «N 

155  XN0RM= XNORM £B (I»J)*B(I»J) 
XN0RM=  SORT (XNORM) 

00  160  I =1 » N 
160  B(I«J)-B(ItJ) /XNORM 
175  CONTINUE 
180  RETURN 
END 
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APPENDIX  B 

SPECIFICATION  OF  PROBLEMS  SOLVED 

Eight  problems  were  specified  by  FluiDyne  for  solution  as 
part  of  the  work  on  this  project.  The  computer  output  for  these 
eight  runs  has  been  delivered  to  FluiDyne.  The  problem  specifi¬ 
cations  for  these  runs  are  given  in  this  appendix.  In  all  solu¬ 
tions  four  divisions  of  the  block  in  the  y  and  z  directions  were 
used.  Ten  division  in  .the  direction  of  the  temperature  gradient 
(the  x  direction)  were  used  in  each  case.  A  value  of  Poisson's 
ratio  -  0.29  was  used  in  each  case. 

Figures  3  and  4  give  respectively  the  variation  of  the 
expansion  coefficient  a  and  the  modulus  E  with  temperature.  The 
variation  of  temperature  with  distance  from  the  heated  face  for 
the  four  conditions  considered  is  shown  in  Fig.  5  through  8. 

The  piecewise  linear  approximations  used  in  the  analyses  are 
also  shown  in  these  figure  .  The  layer  thicknesses  and  proper¬ 
ties  used  in  the  runs  are  given  in  Table  3  through  6. 
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Thermal  Expans 


Fig.  1-4  MODULUS  OF  ELASTICITY,  E  VS  TEMPERATURE  T 
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Temperature 
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Fig. 1-5  PIECEWISE  CONSTANT  TEMPERATURE  APPROXIMATION 
FOR  TIME  «  0.1  SEC  (RUN  1) 
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4000 


3000 


2000 


1000 


Fig.  1-6  PIECEWISE  CONSTANT  TEMPERATURE  APPROXIMATION 
FOR  TIME  *  10  SEC,  THICKNESS  t  =  1  IN. 

(Runs  2,  4,  5  and  7) 
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Distance  From  Heater  Surface,  In. 


Fig.  1-7  PIECEWISE  CONSTANT  TEMPERATURE  APPROXIMATION 

FOR  TIME  »  10  SEC,  THICKNESS  t  -  0.25  IN.  (Run  6) 


Distance  From  Hearer  Surface,  in. 


Fig. 1-8  PIECEWISE  CONSTANT  TEMPERATURE  APPROXIMATION 
FOR  TIME  «?  40  SEC,  (Runs  3  and  8) 
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TABLE  1-2 

PROBLEM  SPECIFICATIONS  FOR  RUNS  1  to  8 


Run 

Condition 

Time 

sec 

Dimensions, 

in. 

a 

b 

t 

1 

Cold  Start 

0.1 

1.0 

.  1.0 

1.0 

2 

Cold  Start 

10.0 

1.0 

1.0 

1.0 

3 

Cold  Start 

40.0 

1.0 

1.0 

1.0 

4 

Cold  Start 

10.0 

2.0 

2.0 

1.0 

5 

Cold  Start 

10.0 

0.5 

0.5 

1.0 

6 

Cold  Start 

10.0 

1.0 

1.0 

0.25 

7 

Cold  Start 

10.0 

1.0 

0.5 

1.0 

8* 

Cold  Start 

40.0 

1.0 

1.0 

1.0 

The  material  properties  for  Run  8  were  constant,  with  the 
value  of  each  constant  at  its  value  at  room  temperature. 
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TABLE  1-3 


DATA  FOR  RUN  1 


Layer 

No. 

Thickness 

H 

Temperature 

T 

Young 1 s 
Modulus 

E 

Thermal 

Coefficient 

a 

1 

0.005 

2480 

9.0xl06 

6. 8x10“ 6 

2 

0.005 

1900 

21.0xl06 

'6. 5x10“  6 

3 

0.005 

1300 

24.5xl06 

6.1x10" 6 

4 

0.010 

650 

26.0xl06 

5.4xl0“6 

5 

0.010 

280 

26.4xl06 

4. 9x10“ 6 

6 

0.015 

140 

26.5xl06 

4. 7x10” 6 

7 

0.025 

80 

26.5xl06 

4. 6x10” 6 

8 

0.125 

60 

26.5xl06 

4. 55x10” 6 

9 

0.300 

60 

26.5xl06 

4. 55x10” 6 

10 

0.500 

60 

26.5xl06 

4. 55x10“ 6 

Cold  Start 
Time  «  0.1  sec 
Thickness  -  1.0  in. 
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TABLE  1-4 

DATA  FOR  RUNS  2,  4,  5,  and  7 


Layer 

No. 

Thickness 

H 

Temperature 

T 

Young's 

Modulus 

E 

Thermal 

Coefficient 

a 

1 

0.020 

3480 

1.8xl06 

7. 2x10" 6 

2 

0.020 

2880 

4.0xl06 

7.0x10" 6 

3 

0.020 

2440 

lO.OxlO6 

6. 8x10" 6 

4 

0.020 

2070 

19.3xl06 

6.6xl0"6 

5 

0.020 

1760 

22 . OxlO6 

6. 4x10" 6 

6 

0.050 

1300 

24.5xl06 

6.1x10" 6 

7 

0.075 

720 

25.9xl06 

5.5xl0"6 

8 

0.150 

260 

26.4xl06 

5. 1x10" 6 

9 

0.225 

80 

26.5xl06 

4.6xl0-6 

10 

0.40 

60 

26.5xl06 

4. 55x10" 6 

Cold  Start 
Time  -  10  sec 
Thickness  »  1.0  in. 
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TABLE  1-5 
DATA  FOR  RUN  6 


Layer 

No. 

Thickness 

H 

Temperature 

T 

Young's 

Modulus 

E 

Thermal 

Coefficient 

a 

1 

0.020 

3480 

1.8xl06 

7.2xl0-6 

2 

0.020 

2880 

4.0xlQ6 

7.0x10” 6 

3 

0.020 

2440 

lO.OxlO6 

6. 8x10" 6 

4 

0.020  . 

2070 

19.3xl06 

6. 6x10” 6 

5 

0.020 

1760 

22.0xl06 

6. 4x10” 6 

6 

0.025 

1440 

23.9xl06 

6. 2x10” 6 

7 

0.025 

1140 

25.0xl06 

5. 9x10” 6 

8 

0.025 

860 

25.6xl06 

5. 6x10“ 6 

9 

0.025 

650 

26.0xl06 

5.4xl0”6 

10 

0.050 

460 

26.2xl06 

5. 2x10" 6 

Cold  Start 
Time  »  10  sec 
Thickness  =  0.25  in. 
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TABLE  1-6 
DATA  FOR  RUN  3 


■  Layer 
No. 

Thickness 

H 

Temperature 

T 

Young ’ s 
Modulus 

E 

Thermal 

Coefficient 

a 

1 

0.025 

3700 

1.4xl06 

7. 2x10" 6 

2 

0.050 

3200 

2.5xl06 

7. 2x10" 6 

3 

0.050 

2600 

6. 8x10 6 

6. 9x10" 6 

4 

0.050 

2100 

19.0xl06 

6. 7x10" 6 

5 

0.075 

1600 

23.2xl06 

6.3xl0"6 

6 

0.100 

1000 

25.4xl06 

5.8xl0"6 

7 

0.150 

540 

26.2xl06 

5 .lxlO"6 

8 

0.150 

220 

26.4xl06 

4.8xl0-6 

9 

0.150 

100 

26.5xl06 

4. 6xl0"6 

10 

0.200 

60 

26.5xl06 

4.55xl0-6 

Cold  Start 
Time  ■  40  sec 
Thickness  ■  1.0  in. 
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APPENDIX  II 

DESCRIPTION  OF  MODIFIED 
THERMAL  STRESS  COMPUTER  PROGRAM,  TPA  II 


Prepared  by: 

Illinois  Institute  of  Technology  Research  Institute 
November  1968 
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DESCRIPTION  OF  TPA  II,  CDC  VERSION 

A.  INTRODUCTION 

TPA  II  is  a  version  of  TPA  which  has  been  modified  to  be 
compatible  with  the  CDC  6600  system  and  to  extend  the  problem 
size  capacity.  The  program  obtains  thermal  stresses  created 
within  a  rectangular  pi....  y  the  variation  of  temperature  and 
material  properties  in  one  direction.  The  governing  equations  are 
formulated  by  means  of  the  finite  element  method  and  then  solved 
by  the  Gauss-Seidel  method.  The  TPA  II  version  can  handle  problems 
with  up  to  20  layers  of  elements  in  the  direction  of  the  thermal 
gradient  and  up  to  nine  elements  in  each  of  the  other  two  directions. 
To  facilitate  convergence  for  these  large  problems,  the  Gauss- 
Seidel  algorithm  has  been  altered  so  as  to  alleviate  the  difficul¬ 
ties  caused  by  the  fixed  points  and  roundoff  error. 

B.  INPUT  DATA 

Let  the  direction  of  the  thermal  gradient  be  denoted  by  x 
and  let  the  remaining  two  orthogonal  directions  be  denoted  by  y 
and  z  respectively.  The  x-direction  must  correspond  to  an  axis 
of  the  prism.  Since  the  prism  is  symmetric  with  respect  to  the 
x-y  and  x-z  planes,  only  a  quadrant  of  the  prism  is  considered. 

The  geometry  of  this  quadrant  is  specified  by  giving  its  y  and  z 
dimensions  and  the  number  of  nodes  desired  in  x,  y,  and  z  direc¬ 
tions.  The  nodes  in  the  y  and  z  directions  will  be  equally  spaced. 
The  spacing  between  the  nodes  in  the  x-direction  which  corresponds 
to  the  thicknesses  of  the  element  layers,  are  specified  by  the 
user  in  the  layer  cards. 
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For  each  layer  of  elements  the  user  must  provide  a  "layer” 
card  which  gives  its  thickness  and  the  temperature,  Young's 
modulus  E,  Poisson's  ratio  v,  and  the  thermal  expansion  coefficient  a. 

The  only  other  required  input  is  composed  of  two  control  num¬ 
bers  NUREAD  and  NUPRNT  .and  the  following  parameters  for  the  Gauss- 
Seidel  method:  the  maximum  number  of  iterations  NIT,  the  error 
printout  cycle  NPERR,  displacement  and  modified  force  printout 
cycle  NPOTPT,  the  tolerance  TOL,  the  initial  relaxation  factor 

i 

XFAC  and  the  error  parameters,  DELERR.  The  initial,  relaxation 
factor  should  be  about  1.5.  After  the  iterative  method  is  started 
if  a  given  iteration  does  not  reduce  the  error  by  a  factor  of 
at  least  DELERR,  the  relaxation  factor  is  reduced.  The  Gauss- 
Seidel  method  terminates  after  the  error  fails  to  decrease, 
when  the  number  of  iterations  exceeds  NIT,  or  if  the  error  meets  the 
tolerance  TOL.  During  the  iterations  the  error  and  displacements 
are  printed  out  according  to  the  value  of  NPERR  and  NPOTPT, 
respectively.  It  is  recommended  that  the  displacements  not  be 
printed  out  at  all  during  the  solution  for  larger  problems  since 
this  output  is  very  time-consuming. 

The  control  number  NUREAD  determines  whether  the  starting  values 
of  the  displacements  are  to  be  read  from  an  input  tape  or  to  be 
computed.  The  second  control  number,  NUPRNT  gives  the  user  the 
option  of  saving  the  displacements  on  tape  or  interpolating  to  a 
lattice  with  a  larger  number  of  elements.  These  options  allow 
the  user  to  run  the  problem  to  a  certain  point,  s’&ve  the  displace¬ 
ments  on  a  tape  and  continue  later  from  the  same  point  by  using  these 
displacements  as  starting  values. 
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The  interpolation  routine  inserts  the  displacement  values 
for  a  lattice  with  n^  elements  in  the  z-direction,  ri2  elements  in 
tae  y-direction,  and  m  elements  in  the  x-direction  to  obtain  the 
starting  values  for  2ri^  and  2^  by  m  lattice.  When  this  option 
is  used  the  second  problem  is  then  automatically  run  with  the 
interpolated  displacements  as  starting  values.  A  detailed 
description  of  the  input  data  required  is  given  in  Table  1. 

C.  DESCRIPTION  OF  OUTPUT 

The  initial  output  is  an  echo  of  all  data.  Thd  final  output 
consists  of  the  displacements  and  stresses.  The  stresses  are  com¬ 
puted  and  printed  at  the  centroids  of  all  elements.  The  stresses 
are  printed  out  in  the  following  manner.  For  each  point  there 
are  five  lines.  The  first  line  gives  the  x,  y  and  z  coordinates 
of  the  point.  The  second  line  gives  the  xx,  yy,  zz,  xy,  xz  and  yz 
components  of  the  stresses  and  the  effective  stress  a&.  The  effec¬ 
tive  stress  is  computed  by  the  formula 

°S.“  (axx-S>2  +  <°yy'S)2+  (|J2Z'S>2+  2(V  +  V  +  ffXZ> 

where 


S  =  1/3(<t  +a  +  .a  ). 

'  v  xx  .  yy  zz' 

The  remaining  three  lines  give  the  three  principal  stresses  and  the 
corresponding  principal  directions.  If  the  solutions  to  additional 
problems  are  desired,  the  data  sets  are  repeated.  Any  number  of 
problems  can  be  solved  in  sequence. 
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tabu; 

INPUT  DATA  FORMAT 


Card  I  TITLE  CARD  (12A6) 

Column 

1-72  Any  alphanumeric  information 
Card  2  PARAMETER  CARD  (6I5.4F10.0) 


1  - 

5 

NX 

Number  of  nodes  in  x-direction  ( 

6  - 

10 

NY 

Number  of  nodes  in  y-direction, 

11  - 

15 

NZ 

Number  of  nodes  in  z-dlrection 

16  - 

20 

NIT 

Maximum  number  of  iterations 

21  - 

25 

NPERR 

Error  printout  cycle 

26  - 

30 

NP0TPT 

Displacement  and  modified  force  printout  cycle 

31  - 

40 

DIMY 

y -dimens ion  of  prism  quadrant 

41  - 

50 

DIMZ 

z-dimension  of  prism  quadrant 

51  - 

60 

TOL 

Error  tolerance  for  Gauss-Scidel 

61  - 

65 

XFAC 

Initial  relaxation  factor 

66  - 

70 

DELERR 

Error  parameter f>,  if  the  error  for  iteration 
i  is  and  the  previous  error  the 

relaxation  factor  is  reduced  if  (-ei_1-e^)/e^<  6. 

Card  3 

Column 

LAYER  PROPERTY  CARD  (I5.5E10.0) 

1.- 

5 

Layer  number 

6  - 

15 

H 

Thickness  of  layer 

16  - 

25 

T 

Temperature  of  layer 

26  - 

35 

E 

Young's  modulus  for  layer 

36  - 

45 

PR 

Poisson's  ratio  for  layer 

46  - 

55 

ALFA 

Thermal  expansion  coefficient  for  layer 

Card  4 

Column 

CONTROL  CARD 

1  - 

5 

NUREAD 

Determines  how  Initial  displacements  are  obtained 

-  0  calculated 

-  2  read  from  tape  8 

6  - 

10 

NUPRNT 

- 

Gives  the  following  options 

■  0  none  of  the  options  used 
-  1  displacements  are  saved  on  tape  10 
m  2  displacements  are  used  to-  interpolate 


to  a  larger  lattice 

Covti'  5  INTERPOLATION  CARD  (OPTIONAL)  (315 ,21110 . 0) 

Column 

1-51  -  Code  number;  must  equal  619’ if  interpolation  is 

to  be  performed. 

6-10  XMAX  -  Number  of  nodes  in  x-direction  for  interpolated  lattice 

11  -  15  JMAX  -  Number  of  nodes  in  x-direction  for  interpolated  lattice 

16  -  20  KMAX  -  Number  of  nodes  in  x-direction  for  interpolated  lattice 

21  -  25  NIT  -  Maximum  number  of  iterations  to  be  performed  for 

second  problem 

26  -  30  NPERR  -  Error  printout  cycle  for  second  problem 

31  -  35  XP0TPT  -  Displacement  printout  cycle  for  second  problem 

36  -  45  T0L  -  Error  tolerance  for  second  problem 

46  -  55  XrAC  -  Initial  relaxation  factor  for  second  problem 

(If  Interpolation  option  is  used,  Card  4  must  be  repeated  after  Card  5) 
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ia.  abstract^  study  ve s  made  of  the  feasibility  of  an  uncooled  throat  for  a  large  hyper¬ 
sonic  wind  tunnel  facility  using  currently  available  materials.  Maximum  fullscale 
stagnation  conditions  would  be  2000psi,  44009R,  1500  lb/sec  air  flow,  and  throat  diam¬ 
eter  10.5-inches.  The  basic  throat  concept  was  that  of  a  ceramic  insulation  layer,  com¬ 
posed  of  small  pieces,  that  would  form  a  protective  liner  within  a  metal  structure. 

High  resistance  to  thermal  spalling  was  the  material  character isticsof  greatest  impor¬ 
tance.  Tests  were  made  of  several  zirconia  materials  and  two  zirconium  diboride  compo¬ 
sitions  by  exposing  them  to  hot  air  flow  in  a  sonic  throat  at  maximum  conditions  of  800 
psi  and  35509R.  Behavior  of  the  zirconia  materials  ranged  from  minor  cracking  to  com¬ 
plete  fragmentation.  The  z ir c onium-dibor ide s  did  not  crack  and  were  oxidation  resis¬ 
tant  at  these  conditions.  In  addition,  the  thermal  stress  distribution  was  studied  for 
the  individual  blocks  that  would  form  the  throat  insulation.  For  this  purpose  the 
three-dimensional  stress  distribution  was  calculated  for  mechanically  unrestrained 
blocks  having  one-dimensional  temperature  distributions.  Effects  of  temperature  dis¬ 
tribution,  block  size  and  block  shape  were  determined.  The  computer  program  is  in¬ 
cluded  with  the  report.  It  was  concluded  that  currently  available  materials  are  not 
satisfactory  for  a  throat  that  would  be  used  with  no  cooling.  Some  of  the  materials 
tested  may  be  satisfactory  if  the  thermal  shock  conditions  were  reduced  by  use  of  film 
drooling  (less  than  that  required  for  a  backside  cooled  throat)  and  by  preheating  with  a 
flow  of  air  through  the  throat  during  heater  pressurization.  A  major  problem  in  use  of 
zirconia  would  he  the  attachment  of  the  insulation  layer  to  the  backup  structure.  Use 
of  zirconium-diboride  would  require  a  design  concept  that  would  he  compatible  with  its 
high  thermal  conductivity.  Both  materials  might  he  used  to  advantage  in  a  single 
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